From 7f02bec9b34e65c17303c1970381cece616c5911 Mon Sep 17 00:00:00 2001 From: Matthew Falcone Date: Sun, 8 Feb 2026 13:10:22 +0000 Subject: [PATCH 1/7] Add mapped inlet changes to new branch --- include/bcs/FoamMappedInletBCBase.h | 70 ++++ include/bcs/FoamMassFlowRateMappedInletBC.h | 18 + include/bcs/FoamScalarBulkMappedInletBC.h | 19 + src/bcs/FoamMappedInletBCBase.C | 389 ++++++++++++++++++ src/bcs/FoamMassFlowRateMappedInletBC.C | 73 ++++ src/bcs/FoamScalarBulkMappedInletBC.C | 74 ++++ test/OpenFOAM/foam_modules.mk | 1 + .../modules/mappedInletTestSolver/Make/files | 3 + .../mappedInletTestSolver/Make/options | 20 + .../mappedInletTestSolver.C | 183 ++++++++ .../mappedInletTestSolver.H | 152 +++++++ test/tests/bcs/mapped_inlet/foam/0/T | 56 +++ test/tests/bcs/mapped_inlet/foam/0/U | 31 ++ test/tests/bcs/mapped_inlet/foam/0/rho | 56 +++ test/tests/bcs/mapped_inlet/foam/constant/g | 21 + .../foam/constant/momentumTransport | 18 + .../foam/constant/physicalProperties | 52 +++ .../mapped_inlet/foam/system/blockMeshDict | 85 ++++ .../foam/system/blockMeshDict.face_edge | 85 ++++ .../foam/system/blockMeshDict.orig | 85 ++++ .../bcs/mapped_inlet/foam/system/controlDict | 44 ++ .../mapped_inlet/foam/system/decomposeParDict | 42 ++ .../foam/system/decomposeParDict.face_edge | 42 ++ .../foam/system/decomposeParDict.orig | 42 ++ .../bcs/mapped_inlet/foam/system/fvSchemes | 51 +++ .../bcs/mapped_inlet/foam/system/fvSolution | 61 +++ test/tests/bcs/mapped_inlet/main.i | 76 ++++ test/tests/bcs/mapped_inlet/main_rotated.i | 50 +++ test/tests/bcs/mapped_inlet/test.py | 91 ++++ test/tests/bcs/mapped_inlet/tests | 153 +++++++ 30 files changed, 2143 insertions(+) create mode 100644 include/bcs/FoamMappedInletBCBase.h create mode 100644 include/bcs/FoamMassFlowRateMappedInletBC.h create mode 100644 include/bcs/FoamScalarBulkMappedInletBC.h create mode 100644 src/bcs/FoamMappedInletBCBase.C create mode 100644 src/bcs/FoamMassFlowRateMappedInletBC.C create mode 100644 src/bcs/FoamScalarBulkMappedInletBC.C create mode 100644 test/OpenFOAM/modules/mappedInletTestSolver/Make/files create mode 100644 test/OpenFOAM/modules/mappedInletTestSolver/Make/options create mode 100644 test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.C create mode 100644 test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.H create mode 100644 test/tests/bcs/mapped_inlet/foam/0/T create mode 100644 test/tests/bcs/mapped_inlet/foam/0/U create mode 100644 test/tests/bcs/mapped_inlet/foam/0/rho create mode 100644 test/tests/bcs/mapped_inlet/foam/constant/g create mode 100644 test/tests/bcs/mapped_inlet/foam/constant/momentumTransport create mode 100644 test/tests/bcs/mapped_inlet/foam/constant/physicalProperties create mode 100644 test/tests/bcs/mapped_inlet/foam/system/blockMeshDict create mode 100644 test/tests/bcs/mapped_inlet/foam/system/blockMeshDict.face_edge create mode 100644 test/tests/bcs/mapped_inlet/foam/system/blockMeshDict.orig create mode 100644 test/tests/bcs/mapped_inlet/foam/system/controlDict create mode 100644 test/tests/bcs/mapped_inlet/foam/system/decomposeParDict create mode 100644 test/tests/bcs/mapped_inlet/foam/system/decomposeParDict.face_edge create mode 100644 test/tests/bcs/mapped_inlet/foam/system/decomposeParDict.orig create mode 100644 test/tests/bcs/mapped_inlet/foam/system/fvSchemes create mode 100644 test/tests/bcs/mapped_inlet/foam/system/fvSolution create mode 100644 test/tests/bcs/mapped_inlet/main.i create mode 100644 test/tests/bcs/mapped_inlet/main_rotated.i create mode 100644 test/tests/bcs/mapped_inlet/test.py create mode 100644 test/tests/bcs/mapped_inlet/tests diff --git a/include/bcs/FoamMappedInletBCBase.h b/include/bcs/FoamMappedInletBCBase.h new file mode 100644 index 00000000..29e31355 --- /dev/null +++ b/include/bcs/FoamMappedInletBCBase.h @@ -0,0 +1,70 @@ +#pragma once + +#include "FoamPostprocessorBCBase.h" +#include + +class FoamMappedInletBCBase : public FoamPostprocessorBCBase +{ +public: + static InputParameters validParams(); + + FoamMappedInletBCBase(const InputParameters & params); + + virtual ~FoamMappedInletBCBase() { destroyCommunicator(_foam_comm); } + +protected: + Foam::vector _offset; + + std::map> _send_map; + + std::map> _recv_map; + + Foam::label _foam_comm; + + MPI_Comm _mpi_comm; + + // create send and receive information for mapping + void createPatchProcMap(); + + // get array from mapped plane on the inlet processes + template + Foam::Field getMappedArray(const Foam::word & name); + + // check if bounding box intersects with rank + bool intersectMapPlane(const Foam::fvMesh & mesh, Real cart_bbox[6]); + + // create/assign communicators for the transfers between map and inlet planes + void createMapComm(const Foam::fvMesh & mesh, + Foam::vectorField face_centres, + std::vector & send_process, + std::vector & recv_process); + + // find index of cell containing point or raise error if not found + int findIndex(const Foam::point & location, const MPI_Comm & comm); + + // handle creation of new communicators in parallel or serial + Foam::label + createCommunicator(const Foam::label parent_comm, std::vector procs, MPI_Comm & new_comm) + { + Foam::label foam_comm; + if (Foam::UPstream::parRun()) + { + Foam::labelList foam_procs(procs.begin(), procs.end()); + foam_comm = Foam::UPstream::allocateCommunicator(parent_comm, foam_procs, true); + new_comm = Foam::PstreamGlobals::MPICommunicators_[foam_comm]; + } + else + { + foam_comm = Foam::UPstream::worldComm; + new_comm = MPI_COMM_WORLD; + } + return foam_comm; + } + + // free communicators if parallel run + void destroyCommunicator(Foam::label comm) + { + if (Foam::UPstream::parRun()) + Foam::UPstream::freeCommunicator(comm); + } +}; diff --git a/include/bcs/FoamMassFlowRateMappedInletBC.h b/include/bcs/FoamMassFlowRateMappedInletBC.h new file mode 100644 index 00000000..39b8520f --- /dev/null +++ b/include/bcs/FoamMassFlowRateMappedInletBC.h @@ -0,0 +1,18 @@ +#pragma once +#include "FoamMappedInletBCBase.h" +#include "MooseEnum.h" + +class FoamMassFlowRateMappedInletBC : public FoamMappedInletBCBase +{ +public: + static InputParameters validParams(); + + FoamMassFlowRateMappedInletBC(const InputParameters & params); + + virtual void imposeBoundaryCondition() override; + +protected: + MooseEnum _scale_method; + + Real _scale_factor; +}; diff --git a/include/bcs/FoamScalarBulkMappedInletBC.h b/include/bcs/FoamScalarBulkMappedInletBC.h new file mode 100644 index 00000000..2a68fd0c --- /dev/null +++ b/include/bcs/FoamScalarBulkMappedInletBC.h @@ -0,0 +1,19 @@ +#pragma once +#include "FoamMappedInletBCBase.h" +#include "MooseEnum.h" + +class FoamScalarBulkMappedInletBC : public FoamMappedInletBCBase +{ +public: + static InputParameters validParams(); + + FoamScalarBulkMappedInletBC(const InputParameters & params); + + virtual void imposeBoundaryCondition() override; + +protected: + MooseEnum _scale_method; + + template + T applyScaleMethod(T & var, const Real bulk_ref, const Real bulk); +}; diff --git a/src/bcs/FoamMappedInletBCBase.C b/src/bcs/FoamMappedInletBCBase.C new file mode 100644 index 00000000..9ca25de2 --- /dev/null +++ b/src/bcs/FoamMappedInletBCBase.C @@ -0,0 +1,389 @@ +#include "FoamMappedInletBCBase.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace +{ +// Get the cartesian bounding box of the mapped inlet plane +void +getBBox(const Foam::vectorField points, Real bbox[6]) +{ + bbox[0] = DBL_MAX; + bbox[1] = DBL_MIN; + bbox[2] = DBL_MAX; + bbox[3] = DBL_MIN; + bbox[4] = DBL_MAX; + bbox[5] = DBL_MIN; + for (auto p : points) + { + bbox[0] = std::min(bbox[0], p.x()); + bbox[1] = std::max(bbox[1], p.x()); + bbox[2] = std::min(bbox[2], p.y()); + bbox[3] = std::max(bbox[3], p.y()); + bbox[4] = std::min(bbox[4], p.z()); + bbox[5] = std::max(bbox[5], p.z()); + } + + MPI_Allreduce(MPI_IN_PLACE, &bbox[0], 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &bbox[1], 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &bbox[2], 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &bbox[3], 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &bbox[4], 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &bbox[5], 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); +} +} + +bool +FoamMappedInletBCBase::intersectMapPlane(const Foam::fvMesh & mesh, Real cart_bbox[6]) +{ + auto & vertices = mesh.points(); + for (int i = 0; i < mesh.nCells(); ++i) + { + auto points = mesh.cellPoints(i); + for (auto pointI : points) + { + // check whether cell intersects bbox + const Foam::point & p = vertices[pointI]; + if (p.x() >= cart_bbox[0] && p.x() <= cart_bbox[1] && p.y() >= cart_bbox[2] && + p.y() <= cart_bbox[3] && p.z() >= cart_bbox[4] && p.z() <= cart_bbox[5]) + { + return true; + } + } + + // The bbox could by narrower than the cell, check in each direction whether bbox is within the + // cell bbox + Real cell_bbox[6] = {DBL_MAX, DBL_MIN, DBL_MAX, DBL_MIN, DBL_MAX, DBL_MIN}; + for (auto point : points) + { + const Foam::point & p = vertices[point]; + cell_bbox[0] = std::min(cell_bbox[0], p.x()); + cell_bbox[1] = std::max(cell_bbox[1], p.x()); + cell_bbox[2] = std::min(cell_bbox[2], p.y()); + cell_bbox[3] = std::max(cell_bbox[3], p.y()); + cell_bbox[4] = std::min(cell_bbox[4], p.z()); + cell_bbox[5] = std::max(cell_bbox[5], p.z()); + } + // check is cart bbox is narrower than cell bbox in x direction + if ((cart_bbox[0] >= cell_bbox[0] && cart_bbox[1] <= cell_bbox[1]) && + (cart_bbox[3] >= cell_bbox[2] || cart_bbox[2] <= cell_bbox[3]) && + (cart_bbox[5] >= cell_bbox[4] || cart_bbox[4] <= cell_bbox[5])) + { + return true; + } + + // check is cart bbox is narrower than cell bbox in y direction + if ((cart_bbox[2] >= cell_bbox[2] && cart_bbox[3] <= cell_bbox[3]) && + (cart_bbox[1] >= cell_bbox[0] || cart_bbox[0] <= cell_bbox[1]) && + (cart_bbox[5] >= cell_bbox[4] || cart_bbox[4] <= cell_bbox[5])) + { + return true; + } + // check is cart bbox is narrower than cell bbox in z direction + if ((cart_bbox[4] >= cell_bbox[4] && cart_bbox[5] <= cell_bbox[5]) && + (cart_bbox[1] >= cell_bbox[0] || cart_bbox[0] <= cell_bbox[1]) && + (cart_bbox[3] >= cell_bbox[2] || cart_bbox[2] <= cell_bbox[3])) + { + return true; + } + } + + return false; +} + +void +FoamMappedInletBCBase::createMapComm(const Foam::fvMesh & mesh, + const Foam::vectorField face_centres, + std::vector & map_process, + std::vector & inlet_process) +{ + Real cart_bbox[6]; + auto mapped_plane = face_centres + _offset; + getBBox(mapped_plane(), cart_bbox); + + int mappedPlaneProcess = intersectMapPlane(mesh, cart_bbox); + int inletPlaneProcess = face_centres.size() > 0; + + std::vector inlet_procs(Foam::UPstream::nProcs()); + std::vector map_procs(Foam::UPstream::nProcs()); + + MPI_Allgather(&mappedPlaneProcess, 1, MPI_INT, map_procs.data(), 1, MPI_INT, MPI_COMM_WORLD); + MPI_Allgather(&inletPlaneProcess, 1, MPI_INT, inlet_procs.data(), 1, MPI_INT, MPI_COMM_WORLD); + + map_process.clear(); + inlet_process.clear(); + + // create list of processes in new communicator and whether they are in the inlet or mapped plane + std::vector processes; + int j = 0; + for (int i = 0; i < Foam::UPstream::nProcs(); ++i) + { + if (inlet_procs[i] || map_procs[i]) + { + processes.push_back(i); + if (inlet_procs[i]) + inlet_process.push_back(j); + if (map_procs[i]) + map_process.push_back(j); + ++j; + } + } + + _foam_comm = createCommunicator(Foam::UPstream::worldComm, processes, _mpi_comm); +} + +void +FoamMappedInletBCBase::createPatchProcMap() +{ + auto & foam_mesh = _mesh->fvMesh(); + auto & boundary = foam_mesh.boundary()[_boundary[0]]; + auto face_centres = boundary.Cf(); + + std::vector map_procs, inlet_procs; + createMapComm(foam_mesh, face_centres, map_procs, inlet_procs); + + if (_mpi_comm == MPI_COMM_NULL) // process not in mapped or inlet planes + return; + + Foam::PstreamBuffers send_points( + Foam::UPstream::commsTypes::nonBlocking, Foam::UPstream::msgType(), _foam_comm); + int rank = Foam::UPstream::myProcNo(_foam_comm); + bool isMapProc = std::find(map_procs.begin(), map_procs.end(), rank) != map_procs.end(); + bool isInletProc = std::find(inlet_procs.begin(), inlet_procs.end(), rank) != inlet_procs.end(); + + if (isInletProc) // send points from inlet process to all map processes + { + for (int proc : map_procs) + { + Foam::UOPstream send(proc, send_points); + send << face_centres; + } + } + + send_points.finishedSends(true); + + std::vector size_requests(inlet_procs.size()); + std::vector data_requests(inlet_procs.size()); + std::vector> recv_indices_procs(inlet_procs.size()); + MPI_Comm map_comm; + auto foam_map_comm = createCommunicator(_foam_comm, map_procs, map_comm); + + if (isMapProc) // check points from each process to see if they are local + { + for (auto i = 0lu; i < inlet_procs.size(); ++i) + { + Foam::vectorField field; + Foam::UIPstream recieve(inlet_procs[i], send_points); + recieve >> field; + auto & vec = _send_map[inlet_procs[i]]; + auto & recv_indices = recv_indices_procs[i]; + for (int j = 0; j < field.size(); ++j) + { + auto index = findIndex(field[j] + _offset, map_comm); + if (index >= 0) + { + vec.push_back(index); // assign to send map required indices + recv_indices.push_back(j); + } + } + if (vec.size() == 0) + _send_map.erase(inlet_procs[i]); + + // Let original processes know which points will come from each rank + int size = recv_indices.size(); + MPI_Isend(&size, 1, MPI_INT, inlet_procs[i], 0, _mpi_comm, &size_requests.at(i)); + MPI_Isend(recv_indices.data(), + recv_indices.size(), + MPI_INT, + inlet_procs[i], + 1, + _mpi_comm, + &data_requests.at(i)); + } + } + + destroyCommunicator(foam_map_comm); + + if (isInletProc) // create map to determine where data from map processes should go + { + for (auto & proc : map_procs) + { + int size; + MPI_Recv(&size, 1, MPI_INT, proc, 0, _mpi_comm, MPI_STATUS_IGNORE); + + std::vector recv_indices(size); + MPI_Recv(recv_indices.data(), size, MPI_INT, proc, 1, _mpi_comm, MPI_STATUS_IGNORE); + for (auto & index : recv_indices) + { + assert(index < face_centres.size()); + _recv_map[proc].push_back(index); + } + } + } + MPI_Barrier(_mpi_comm); +} + +int +FoamMappedInletBCBase::findIndex(const Foam::point & location, const MPI_Comm & comm) +{ + /* + This function uses several ways of finding the mapped plane cell + 1. use findCell function + - Sometimes on cell boundaries it may not find the cell or two processes with both find it + 2. If none found, find the closest point and do an expanded bounding box search, raise error + if none still found. + 3. If multiple found, use the cell with cell centre closer to desired location + 4. If still multiple found, use the cell closer to the inlet. + - In the unlikely chance there are still multiple cells detected, raise a warning + Note that it is possible there is some non-deterministic behaviour in the function but this + shouldn't be a problem in practice. + */ + int index = _mesh->fvMesh().findCell(location, Foam::polyMesh::FACE_PLANES); + int gl_index; + MPI_Allreduce(&index, &gl_index, 1, MPI_INT, MPI_MAX, comm); + + // expand cell bounding box and repeat search + if (gl_index < 0) + { + Foam::label celli = _mesh->fvMesh().findNearestCell(location); + + bool in_cell = _mesh->fvMesh().pointInCellBB(location, celli, 0.1); + index = (in_cell) ? celli : -1; + MPI_Allreduce(&index, &gl_index, 1, MPI_INT, MPI_MAX, comm); + + int rank; + MPI_Comm_rank(comm, &rank); + if (gl_index < 0 && rank == 0) + { + mooseError("Face centre at location (", + location[0], + ",", + location[1], + ",", + location[2], + ") does not have a mapped plane location"); + } + } + + // use cell with cell centre closest to the location + Foam::scalar dist{DBL_MAX}, gl_dist; + if (index != -1) + dist = Foam::mag(_mesh->fvMesh().cellCentres()[index] - location); + + MPI_Allreduce(&dist, &gl_dist, 1, MPI_DOUBLE, MPI_MIN, comm); + if (dist != gl_dist) + index = -1; + + // 2. use cell centre closest to inlet point + int in_cell = index != -1; + MPI_Allreduce(MPI_IN_PLACE, &in_cell, 1, MPI_INT, MPI_SUM, comm); + if (in_cell > 1) + { + if (index != -1) + dist = Foam::mag(_mesh->fvMesh().cellCentres()[index] - (location - _offset)); + MPI_Allreduce(&dist, &gl_dist, 1, MPI_DOUBLE, MPI_MIN, comm); + if (dist != gl_dist) + index = -1; + } + in_cell = index != -1; + MPI_Allreduce(MPI_IN_PLACE, &in_cell, 1, MPI_INT, MPI_SUM, comm); + + if (in_cell > 1) + mooseWarning("More than 1 process found location (", + location[0], + ",", + location[1], + ",", + location[2], + ")"); + return index; +} + +InputParameters +FoamMappedInletBCBase::validParams() +{ + auto params = FoamPostprocessorBCBase::validParams(); + params.addRequiredParam>("translation_vector", + "A vector indicating the location of recycling plane"); + + return params; +} + +FoamMappedInletBCBase::FoamMappedInletBCBase(const InputParameters & params) + : FoamPostprocessorBCBase(params), + _offset(), + _send_map(), + _recv_map(), + _foam_comm(0), + _mpi_comm(MPI_COMM_NULL) +{ + if (_boundary.size() > 1) + mooseError("There can only be one boundary using this method"); + + auto param_offset = params.get>("translation_vector"); + assert(param_offset.size() == 3); + + _offset = {param_offset[0], param_offset[1], param_offset[2]}; + createPatchProcMap(); +} + +template +Foam::Field +FoamMappedInletBCBase::getMappedArray(const Foam::word & name) +{ + if (_mpi_comm == MPI_COMM_NULL) + return Foam::Field(); + + auto & foam_mesh = _mesh->fvMesh(); + auto & boundary_patch = foam_mesh.boundary()[_boundary[0]]; + + Foam::PstreamBuffers sendBuf( + Foam::UPstream::commsTypes::nonBlocking, Foam::UPstream::msgType(), _foam_comm); + if (_send_map.size() > 0) + { + auto & var = foam_mesh.lookupObject>(name); + for (auto & pair : _send_map) + { + auto & proc = pair.first; + auto & send_indices = pair.second; + + Foam::Field points(send_indices.size()); + for (auto j = 0lu; j < send_indices.size(); ++j) + points[j] = var[send_indices[j]]; + + Foam::UOPstream send(proc, sendBuf); + send << points; + } + } + sendBuf.finishedSends(true); + + Foam::Field boundaryData(boundary_patch.size()); + if (_recv_map.size() > 0) + { + for (auto & pair : _recv_map) + { + auto & proc = pair.first; + auto & recv_indices = pair.second; + + Foam::UIPstream recv(proc, sendBuf); + Foam::Field recvData; + recv >> recvData; + for (auto j = 0lu; j < recv_indices.size(); ++j) + { + boundaryData[recv_indices[j]] = recvData[j]; + } + } + } + + return boundaryData; +} diff --git a/src/bcs/FoamMassFlowRateMappedInletBC.C b/src/bcs/FoamMassFlowRateMappedInletBC.C new file mode 100644 index 00000000..88ba1f96 --- /dev/null +++ b/src/bcs/FoamMassFlowRateMappedInletBC.C @@ -0,0 +1,73 @@ +#include "FoamMassFlowRateMappedInletBC.h" +#include "InputParameters.h" +#include "MooseTypes.h" +#include "PstreamReduceOps.H" +#include "Registry.h" + +#include "ops.H" +#include "vectorField.H" +#include "volFieldsFwd.H" + +registerMooseObject("hippoApp", FoamMassFlowRateMappedInletBC); + +InputParameters +FoamMassFlowRateMappedInletBC::validParams() +{ + auto params = FoamMappedInletBCBase::validParams(); + + params.addParam("scale_factor", 1., "Scale factor multiply mass flow rate pp by."); + MooseEnum scaleEnum("SCALE NONE", "SCALE"); + params.addParam("scale_method", + scaleEnum, + "Method used to maintain inlet bulk variable. " + "SCALE means the variable is multiplied by a factor, " + "NONE means the variable is not scaled."); + + params.remove("foam_variable"); + params.addPrivateParam("foam_variable", "U"); + + return params; +} + +FoamMassFlowRateMappedInletBC::FoamMassFlowRateMappedInletBC(const InputParameters & params) + : FoamMappedInletBCBase(params), + _scale_method(params.get("scale_method")), + _scale_factor(params.get("scale_factor")) +{ +} + +void +FoamMassFlowRateMappedInletBC::imposeBoundaryCondition() +{ + auto & foam_mesh = _mesh->fvMesh(); + auto & boundary_patch = foam_mesh.boundary()[_boundary[0]]; + + // currently we map mass flux rather than velocity, maybe useful to have option + auto && U_map = getMappedArray("U"); + auto && rho_map = getMappedArray("rho"); + auto g_map = rho_map * U_map; + + auto & rho = boundary_patch.lookupPatchField("rho"); + auto & U_var = const_cast &>( + boundary_patch.lookupPatchField("U")); + + Foam::vectorField g_var(U_var.size()); + g_var = rho_map * U_map; + if (_scale_method == "SCALE") + { + auto & Sf = boundary_patch.Sf(); + + auto m_dot = Foam::returnReduce(Foam::sum(g_map & -Sf), Foam::sumOp()); + if (fabs(m_dot) > 1e-8) + { + g_var *= _scale_factor * _pp_value / m_dot; + } + else + { + auto area = Foam::returnReduce(Foam::sum(boundary_patch.magSf()), Foam::sumOp()); + g_var -= _scale_factor * boundary_patch.nf() * _pp_value / area; + } + } + + U_var == g_var / rho; +} diff --git a/src/bcs/FoamScalarBulkMappedInletBC.C b/src/bcs/FoamScalarBulkMappedInletBC.C new file mode 100644 index 00000000..f1cb83d9 --- /dev/null +++ b/src/bcs/FoamScalarBulkMappedInletBC.C @@ -0,0 +1,74 @@ +#include "FoamScalarBulkMappedInletBC.h" +#include "InputParameters.h" +#include "MooseTypes.h" +#include "PstreamReduceOps.H" +#include "Registry.h" + +#include "ops.H" +#include "volFieldsFwd.H" + +registerMooseObject("hippoApp", FoamScalarBulkMappedInletBC); + +InputParameters +FoamScalarBulkMappedInletBC::validParams() +{ + auto params = FoamMappedInletBCBase::validParams(); + MooseEnum scaleEnum("SCALE SUBTRACT NONE", "SCALE"); + params.addParam("scale_method", + scaleEnum, + "Method used to maintain inlet bulk variable. " + "SCALE means the variable is multiplied by a factor, " + "SUBTRACT means the variable is reduced by constant," + "NONE means the variable is not scaled."); + + return params; +} + +FoamScalarBulkMappedInletBC::FoamScalarBulkMappedInletBC(const InputParameters & params) + : FoamMappedInletBCBase(params), _scale_method(params.get("scale_method")) +{ +} + +void +FoamScalarBulkMappedInletBC::imposeBoundaryCondition() +{ + auto & foam_mesh = _mesh->fvMesh(); + auto & boundary_patch = foam_mesh.boundary()[_boundary[0]]; + + auto && var_map = getMappedArray(_foam_variable); + auto & Sf = boundary_patch.magSf(); + + auto totalArea = Foam::sum(Sf); + Foam::reduce(totalArea, Foam::sumOp()); + + auto var_bulk = Foam::sum(var_map * Sf) / totalArea; + + Foam::reduce(var_bulk, Foam::sumOp()); + + auto & var = const_cast &>( + boundary_patch.lookupPatchField(_foam_variable)); + + var == applyScaleMethod(var_map, _pp_value, var_bulk); +} + +template +T +FoamScalarBulkMappedInletBC::applyScaleMethod(T & var, const Real bulk_ref, const Real bulk) +{ + if (_scale_method == "SCALE") + { + return (var * bulk_ref / bulk)(); + } + else if (_scale_method == "SUBTRACT") + { + return (var + bulk_ref - bulk)(); + } + else if (_scale_method == "NONE") + { + return var; + } + else + { + mooseError("Invalid scale method '", _scale_method, "'."); + } +} diff --git a/test/OpenFOAM/foam_modules.mk b/test/OpenFOAM/foam_modules.mk index 10eef08a..6ce6c3fc 100644 --- a/test/OpenFOAM/foam_modules.mk +++ b/test/OpenFOAM/foam_modules.mk @@ -13,3 +13,4 @@ build_foam_tests: +@$(WMAKE) -s -j $(MOOSE_JOBS) test/OpenFOAM/modules/laplacianTestSolver/ +@$(WMAKE) -s -j $(MOOSE_JOBS) test/OpenFOAM/modules/odeTestSolver/ +@$(WMAKE) -s -j $(MOOSE_JOBS) test/OpenFOAM/modules/postprocessorTestSolver/ + +@$(WMAKE) -s -j $(MOOSE_JOBS) test/OpenFOAM/modules/mappedInletTestSolver/ diff --git a/test/OpenFOAM/modules/mappedInletTestSolver/Make/files b/test/OpenFOAM/modules/mappedInletTestSolver/Make/files new file mode 100644 index 00000000..0f5991d4 --- /dev/null +++ b/test/OpenFOAM/modules/mappedInletTestSolver/Make/files @@ -0,0 +1,3 @@ +SOURCE += mappedInletTestSolver.C + +LIB = $(FOAM_USER_LIBBIN)/libmappedInletTestSolver diff --git a/test/OpenFOAM/modules/mappedInletTestSolver/Make/options b/test/OpenFOAM/modules/mappedInletTestSolver/Make/options new file mode 100644 index 00000000..cc845ca8 --- /dev/null +++ b/test/OpenFOAM/modules/mappedInletTestSolver/Make/options @@ -0,0 +1,20 @@ +EXE_INC = \ + -I$(LIB_SRC)/physicalProperties/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \ + -I$(LIB_SRC)/ThermophysicalTransportModels/thermophysicalTransportModel/lnInclude \ + -I$(LIB_SRC)/ThermophysicalTransportModels/solid/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + +LIB_LIBS = \ + -lsolidThermo \ + -lsolidThermophysicalTransportModels \ + -lcoupledThermophysicalTransportModels \ + -lspecie \ + -lfiniteVolume \ + -lmeshTools \ + -lsampling \ + -lfvModels \ + -lfvConstraints diff --git a/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.C b/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.C new file mode 100644 index 00000000..17f910c8 --- /dev/null +++ b/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.C @@ -0,0 +1,183 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2022-2024 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "dimensionSets.H" +#include "dimensionedScalar.H" +#include "dimensionedVector.H" +#include "fvMesh.H" +#include "mappedInletTestSolver.H" +#include "fvMeshMover.H" +#include "addToRunTimeSelectionTable.H" +#include "fvConstraints.H" +#include "fvmLaplacian.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solvers +{ +defineTypeNameAndDebug(mappedInletTestSolver, 0); +addToRunTimeSelectionTable(solver, mappedInletTestSolver, fvMesh); +} +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // + +bool +Foam::solvers::mappedInletTestSolver::dependenciesModified() const +{ + return runTime.controlDict().modified(); +} + +bool +Foam::solvers::mappedInletTestSolver::read() +{ + solver::read(); + + maxDeltaT_ = runTime.controlDict().found("maxDeltaT") + ? runTime.controlDict().lookup("maxDeltaT", runTime.userUnits()) + : vGreat; + + return true; +} + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // +// Solver based on solid.C module +Foam::solvers::mappedInletTestSolver::mappedInletTestSolver(fvMesh & mesh, + autoPtr thermoPtr) + : solver(mesh), + + thermoPtr_(thermoPtr), + thermo_(thermoPtr_()), + + T_(IOobject("T", mesh.time().name(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh), + + U_(IOobject("U", mesh.time().name(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh), + + thermophysicalTransport(solidThermophysicalTransportModel::New(thermo_)), + thermo(thermo_), + T(T_), + U(U_) +{ + thermo.validate("solid", "h", "e"); +} + +Foam::solvers::mappedInletTestSolver::mappedInletTestSolver(fvMesh & mesh) + : mappedInletTestSolver(mesh, solidThermo::New(mesh)) +{ + // Read the controls + read(); +} + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::solvers::mappedInletTestSolver::~mappedInletTestSolver() {} + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::scalar +Foam::solvers::mappedInletTestSolver::maxDeltaT() const +{ + return min(fvModels().maxDeltaT(), maxDeltaT_); +} + +void +Foam::solvers::mappedInletTestSolver::preSolve() +{ + fvModels().preUpdateMesh(); + + // Update the mesh for topology change, mesh to mesh mapping + mesh_.update(); +} + +void +Foam::solvers::mappedInletTestSolver::moveMesh() +{ + if (pimple.firstIter() || pimple.moveMeshOuterCorrectors()) + { + if (!mesh_.mover().solidBody()) + { + FatalErrorInFunction << "Region " << name() << " of type " << type() + << " does not support non-solid body mesh motion" << exit(FatalError); + } + + mesh_.move(); + } +} + +void +Foam::solvers::mappedInletTestSolver::motionCorrector() +{ +} + +void +Foam::solvers::mappedInletTestSolver::prePredictor() +{ +} + +void +Foam::solvers::mappedInletTestSolver::momentumPredictor() +{ + auto & coords = mesh.C().primitiveField(); + auto & U_field = U_.primitiveFieldRef(); + auto time = mesh.time().userTimeValue(); + + auto x = coords.component(0)(); + auto y = coords.component(1)(); + auto z = coords.component(2)(); + + U_field.replace(0, x + y + z + time); + U_field.replace(1, x - y + z + time); + U_field.replace(2, x + y - z + time); +} + +void +Foam::solvers::mappedInletTestSolver::thermophysicalPredictor() +{ + volScalarField & e = thermo.he(); + auto & coords = mesh.C().primitiveField(); + e.primitiveFieldRef() = mag(coords) + mesh.time().userTimeValue(); + thermo_.correct(); +} + +void +Foam::solvers::mappedInletTestSolver::pressureCorrector() +{ +} + +void +Foam::solvers::mappedInletTestSolver::postCorrector() +{ +} + +void +Foam::solvers::mappedInletTestSolver::postSolve() +{ +} + +// ************************************************************************* // diff --git a/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.H b/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.H new file mode 100644 index 00000000..24705868 --- /dev/null +++ b/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.H @@ -0,0 +1,152 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2022-2023 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::solvers::mappedInletTestSolver + +Description + Solver module for to test mapped inlet implementation method + +SourceFiles + mappedInletTestSolver.C + +\*---------------------------------------------------------------------------*/ + +#pragma once + +#include "solver.H" +#include "solidThermophysicalTransportModel.H" +#include "volFieldsFwd.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solvers +{ + +/*---------------------------------------------------------------------------*\ + Class mappedInletTestSolver Declaration +\*---------------------------------------------------------------------------*/ + +class mappedInletTestSolver : public solver +{ + +protected: + // Control parameters + scalar maxDeltaT_; + + // Thermophysical properties + autoPtr thermoPtr_; + + solidThermo & thermo_; + + volScalarField T_; + + volVectorField U_; + + autoPtr thermophysicalTransport; + + // Protected Member Functions + + //- Return true if the solver's dependencies have been modified + virtual bool dependenciesModified() const; + + //- Read controls + virtual bool read(); + +public: + // Public Data + // reference to thermophysical properties + solidThermo & thermo; + + //- Reference to the temperature field + const volScalarField & T; + + //- Reference velocity field + const volVectorField & U; + + //- Runtime type information + TypeName("mappedInletTestSolver"); + + // Constructors + + //- Construct from region mesh + mappedInletTestSolver(fvMesh & mesh, autoPtr thermoPtr); + + mappedInletTestSolver(fvMesh & mesh); + + //- Disallow default bitwise copy construction + mappedInletTestSolver(const mappedInletTestSolver &) = delete; + + //- Destructor + virtual ~mappedInletTestSolver(); + + // Member Functions + + //- Return the current maximum time-step for stable solution + virtual scalar maxDeltaT() const; + + //- Called at the start of the time-step, before the PIMPLE loop + virtual void preSolve(); + + //- Called at the start of the PIMPLE loop to move the mesh + virtual void moveMesh(); + + //- Corrections that follow mesh motion + virtual void motionCorrector(); + + //- Called at the beginning of the PIMPLE loop + virtual void prePredictor(); + + //- Construct and optionally solve the momentum equation + virtual void momentumPredictor(); + + //- Construct and solve the energy equation, + // convert to temperature + // and update thermophysical and transport properties + virtual void thermophysicalPredictor(); + + //- Construct and solve the pressure equation in the PISO loop + virtual void pressureCorrector(); + + //- Correct the thermophysical transport modelling + virtual void postCorrector(); + + //- Called after the PIMPLE loop at the end of the time-step + virtual void postSolve(); + + // Member Operators + + //- Disallow default bitwise assignment + void operator=(const mappedInletTestSolver &) = delete; +}; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solvers +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/foam/0/T b/test/tests/bcs/mapped_inlet/foam/0/T new file mode 100644 index 00000000..69c3f725 --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/0/T @@ -0,0 +1,56 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 10 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class volScalarField; + location "0"; + object T; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 1 0 0 0]; + +internalField uniform 2.; + +boundaryField +{ + left + { + type fixedValue; + value uniform 0.; + } + right + { + type calculated; + value uniform 0.; + } + top + { + type calculated; + value uniform 0.; + } + bottom + { + type fixedValue; + value uniform 0.; + } + back + { + type calculated; + value uniform 0.; + } + front + { + type fixedValue; + value uniform 0.; + } +} + + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/foam/0/U b/test/tests/bcs/mapped_inlet/foam/0/U new file mode 100644 index 00000000..15b296fa --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/0/U @@ -0,0 +1,31 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 12 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class volVectorField; + location "0"; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [ 0 1 -1 0 0 0 0 ]; + +internalField uniform (1 -0.5 0.25); + +boundaryField +{ + + ".*" + { + type fixedValue; + value $internalField; + } +} + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/foam/0/rho b/test/tests/bcs/mapped_inlet/foam/0/rho new file mode 100644 index 00000000..d30a04dc --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/0/rho @@ -0,0 +1,56 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 12 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class volScalarField; + location "0"; + object rho; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -3 0 0 0 0 0]; + +internalField uniform 1; + +boundaryField +{ + left + { + type calculated; + value uniform 1; + } + right + { + type calculated; + value uniform 1; + } + top + { + type calculated; + value uniform 1; + } + front + { + type calculated; + value uniform 1; + } + bottom + { + type calculated; + value uniform 1; + } + back + { + type calculated; + value uniform 1; + } +} + + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/foam/constant/g b/test/tests/bcs/mapped_inlet/foam/constant/g new file mode 100644 index 00000000..8af96f3a --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/constant/g @@ -0,0 +1,21 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 10 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class uniformDimensionedVectorField; + location "constant"; + object g; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -2 0 0 0 0]; +value (0 0 0); + + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/foam/constant/momentumTransport b/test/tests/bcs/mapped_inlet/foam/constant/momentumTransport new file mode 100644 index 00000000..0416f1a9 --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/constant/momentumTransport @@ -0,0 +1,18 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 10 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + object RASProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType laminar; + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/foam/constant/physicalProperties b/test/tests/bcs/mapped_inlet/foam/constant/physicalProperties new file mode 100644 index 00000000..4c5dac6a --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/constant/physicalProperties @@ -0,0 +1,52 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 10 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "constant"; + object physicalProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType +{ + type heSolidThermo; + mixture pureMixture; + transport constIsoSolid; + thermo eConst; + equationOfState rhoConst; + specie specie; + energy sensibleInternalEnergy; +} + +mixture +{ + specie + { + molWeight 1; + } + thermodynamics + { + Cv 1; // Specific heat capacity [J/(kg·K)] + Hf 1; // Heat of formation [J/kg] + Tref 0; + } + transport + { + kappa 1; // Thermal conductivity [W/(m·K)] + mu 1.; + } + equationOfState + { + rho 0.5; // Density [kg/m^3] + } +} + + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/foam/system/blockMeshDict b/test/tests/bcs/mapped_inlet/foam/system/blockMeshDict new file mode 100644 index 00000000..8221a4b1 --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/system/blockMeshDict @@ -0,0 +1,85 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} + +vertices +( + ( 0.0 0.0 0.0 ) + ( 2.0 0.0 0.0 ) + ( 2.0 2.0 0.0 ) + ( 0.0 2.0 0.0 ) + + ( 0.0 0.0 2.0) + ( 2.0 0.0 2.0) + ( 2.0 2.0 2.0) + ( 0.0 2.0 2.0) +); + +blocks +( + hex (0 1 2 3 4 5 6 7) (9 9 9) simpleGrading (1 1 1) +); + +boundary +( + + // interface + left + { + type wall; + faces + ( + (4 7 3 0) + ); + } + + right + { + type wall; + faces + ( + (6 5 1 2) + ); + } + + top + { + type wall; + faces + ( + (2 3 7 6) + ); + } + + front + { + type wall; + faces + ( + (3 2 1 0) + ); + } + + bottom + { + type wall; + faces + ( + (0 1 5 4) + ); + } + + back + { + type wall; + faces + ( + (4 5 6 7) + ); + } + +); diff --git a/test/tests/bcs/mapped_inlet/foam/system/blockMeshDict.face_edge b/test/tests/bcs/mapped_inlet/foam/system/blockMeshDict.face_edge new file mode 100644 index 00000000..4a877e58 --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/system/blockMeshDict.face_edge @@ -0,0 +1,85 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} + +vertices +( + ( 0.0 0.0 0.0 ) + ( 2.0 0.0 0.0 ) + ( 2.0 2.0 0.0 ) + ( 0.0 2.0 0.0 ) + + ( 0.0 0.0 2.0) + ( 2.0 0.0 2.0) + ( 2.0 2.0 2.0) + ( 0.0 2.0 2.0) +); + +blocks +( + hex (0 1 2 3 4 5 6 7) (8 8 8) simpleGrading (1 1 1) +); + +boundary +( + + // interface + left + { + type wall; + faces + ( + (4 7 3 0) + ); + } + + right + { + type wall; + faces + ( + (6 5 1 2) + ); + } + + top + { + type wall; + faces + ( + (2 3 7 6) + ); + } + + front + { + type wall; + faces + ( + (3 2 1 0) + ); + } + + bottom + { + type wall; + faces + ( + (0 1 5 4) + ); + } + + back + { + type wall; + faces + ( + (4 5 6 7) + ); + } + +); diff --git a/test/tests/bcs/mapped_inlet/foam/system/blockMeshDict.orig b/test/tests/bcs/mapped_inlet/foam/system/blockMeshDict.orig new file mode 100644 index 00000000..8221a4b1 --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/system/blockMeshDict.orig @@ -0,0 +1,85 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} + +vertices +( + ( 0.0 0.0 0.0 ) + ( 2.0 0.0 0.0 ) + ( 2.0 2.0 0.0 ) + ( 0.0 2.0 0.0 ) + + ( 0.0 0.0 2.0) + ( 2.0 0.0 2.0) + ( 2.0 2.0 2.0) + ( 0.0 2.0 2.0) +); + +blocks +( + hex (0 1 2 3 4 5 6 7) (9 9 9) simpleGrading (1 1 1) +); + +boundary +( + + // interface + left + { + type wall; + faces + ( + (4 7 3 0) + ); + } + + right + { + type wall; + faces + ( + (6 5 1 2) + ); + } + + top + { + type wall; + faces + ( + (2 3 7 6) + ); + } + + front + { + type wall; + faces + ( + (3 2 1 0) + ); + } + + bottom + { + type wall; + faces + ( + (0 1 5 4) + ); + } + + back + { + type wall; + faces + ( + (4 5 6 7) + ); + } + +); diff --git a/test/tests/bcs/mapped_inlet/foam/system/controlDict b/test/tests/bcs/mapped_inlet/foam/system/controlDict new file mode 100644 index 00000000..2bbe0e5f --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/system/controlDict @@ -0,0 +1,44 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 10 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solver mappedInletTestSolver; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 5.; + +deltaT 1.; + +writeControl timeStep; + +writeInterval 1; + +writeFormat ascii; + +writePrecision 20; + +writeCompression off; + +timeFormat general; + +timePrecision 20; + +runTimeModifiable true; + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/foam/system/decomposeParDict b/test/tests/bcs/mapped_inlet/foam/system/decomposeParDict new file mode 100644 index 00000000..9df2378c --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/system/decomposeParDict @@ -0,0 +1,42 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 10 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "system"; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 5; + +method scotch; + +simpleCoeffs +{ + n (2 1 1); +} + +hierarchicalCoeffs +{ + n (1 1 1); + order xyz; +} + +manualCoeffs +{ + dataFile ""; +} + +distributed no; + +roots ( ); + + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/foam/system/decomposeParDict.face_edge b/test/tests/bcs/mapped_inlet/foam/system/decomposeParDict.face_edge new file mode 100644 index 00000000..0204a6b9 --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/system/decomposeParDict.face_edge @@ -0,0 +1,42 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 10 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "system"; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 2; + +method simple; + +simpleCoeffs +{ + n (2 1 1); +} + +hierarchicalCoeffs +{ + n (1 1 1); + order xyz; +} + +manualCoeffs +{ + dataFile ""; +} + +distributed no; + +roots ( ); + + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/foam/system/decomposeParDict.orig b/test/tests/bcs/mapped_inlet/foam/system/decomposeParDict.orig new file mode 100644 index 00000000..9df2378c --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/system/decomposeParDict.orig @@ -0,0 +1,42 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 10 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "system"; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 5; + +method scotch; + +simpleCoeffs +{ + n (2 1 1); +} + +hierarchicalCoeffs +{ + n (1 1 1); + order xyz; +} + +manualCoeffs +{ + dataFile ""; +} + +distributed no; + +roots ( ); + + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/foam/system/fvSchemes b/test/tests/bcs/mapped_inlet/foam/system/fvSchemes new file mode 100644 index 00000000..787f3b83 --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/system/fvSchemes @@ -0,0 +1,51 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 12 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss linear; +} + +divSchemes +{ + default none; + + div(phi,U) Gauss linear; + div(phi,K) Gauss linear; + div(phi,h) Gauss linear; + div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; +} + +laplacianSchemes +{ + default Gauss linear corrected; +} + +interpolationSchemes +{ + default linear; +} + +snGradSchemes +{ + default corrected; +} + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/foam/system/fvSolution b/test/tests/bcs/mapped_inlet/foam/system/fvSolution new file mode 100644 index 00000000..61143b21 --- /dev/null +++ b/test/tests/bcs/mapped_inlet/foam/system/fvSolution @@ -0,0 +1,61 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 12 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + rho + { + solver diagonal; + } + + rhoFinal + { + $rho; + } + + + "(U|h|p_rgh)" + { + solver PBiCGStab; + preconditioner DILU; + tolerance 1e-8; + relTol 1e-8; + } + + "(U|h|p_rgh)Final" + { + $U; + tolerance 1e-8; + relTol 1e-8; + } +} + +PIMPLE +{ + momentumPredictor yes; + pRefCell 0; + pRefValue 0; +} + +relaxationFactors +{ + equations + { + h 1; + U 1; + } +} + +// ************************************************************************* // diff --git a/test/tests/bcs/mapped_inlet/main.i b/test/tests/bcs/mapped_inlet/main.i new file mode 100644 index 00000000..50a34648 --- /dev/null +++ b/test/tests/bcs/mapped_inlet/main.i @@ -0,0 +1,76 @@ +[Mesh] + type = FoamMesh + case = 'foam' + foam_patch = 'left right bottom top back front' +[] + +[FoamBCs] + [mass_flowx] + type=FoamMassFlowRateMappedInletBC + boundary = 'left' + default = 1 + translation_vector = '1. 0 0' + [] + [tempx] + type=FoamScalarBulkMappedInletBC + boundary = 'left' + default = 1 + translation_vector = '1. 0 0' + foam_variable = 'T' + [] + [mass_flowy] + type=FoamMassFlowRateMappedInletBC + boundary = 'bottom' + default = 1 + translation_vector = '0 1. 0' + [] + [tempy] + type=FoamScalarBulkMappedInletBC + boundary = 'bottom' + default = 1 + translation_vector = '0 1. 0' + foam_variable = 'T' + [] + [mass_flowz] + type=FoamMassFlowRateMappedInletBC + boundary = 'front' + default = 1 + translation_vector = '0 0 1.' + [] + [tempz] + type=FoamScalarBulkMappedInletBC + boundary = 'front' + default = 1 + translation_vector = '0 0 1.' + foam_variable = 'T' + [] +[] + +[Postprocessors] + [pp] + type = ParsedPostprocessor + expression = '2' + execute_on = TIMESTEP_BEGIN + [] +[] + + +[Problem] + type = FoamProblem + # Take the boundary temperature from OpenFOAM and set it on the MOOSE mesh. +[] + +[Executioner] + type = Transient + end_time = 5 + [TimeSteppers] + [foam] + type = FoamControlledTimeStepper + [] + [] +[] + +[Outputs] + exodus = true + csv=true +[] diff --git a/test/tests/bcs/mapped_inlet/main_rotated.i b/test/tests/bcs/mapped_inlet/main_rotated.i new file mode 100644 index 00000000..553b6d3d --- /dev/null +++ b/test/tests/bcs/mapped_inlet/main_rotated.i @@ -0,0 +1,50 @@ +[Mesh] + type = FoamMesh + case = 'foam' + foam_patch = 'left right bottom top back front' +[] + +[FoamBCs] + [mass_flowx] + type=FoamMassFlowRateMappedInletBC + boundary = 'left' + default = 1 + translation_vector = '${fparse sqrt(0.5)} ${fparse sqrt(0.5)} 0' + [] + [tempx] + type=FoamScalarBulkMappedInletBC + boundary = 'left' + default = 1 + translation_vector = '${fparse sqrt(0.5)} ${fparse sqrt(0.5)} 0' + foam_variable = 'T' + [] +[] + +[Postprocessors] + [pp] + type = ParsedPostprocessor + expression = '2' + execute_on = TIMESTEP_BEGIN + [] +[] + + +[Problem] + type = FoamProblem + # Take the boundary temperature from OpenFOAM and set it on the MOOSE mesh. +[] + +[Executioner] + type = Transient + end_time = 5 + [TimeSteppers] + [foam] + type = FoamControlledTimeStepper + [] + [] +[] + +[Outputs] + exodus = true + csv=true +[] diff --git a/test/tests/bcs/mapped_inlet/test.py b/test/tests/bcs/mapped_inlet/test.py new file mode 100644 index 00000000..015437dc --- /dev/null +++ b/test/tests/bcs/mapped_inlet/test.py @@ -0,0 +1,91 @@ +"""Tests for imposing BCs in OpenFOAM using MOOSE input file syntax""" + +import unittest +import fluidfoam as ff +import numpy as np + +from read_hippo_data import get_foam_times + +CASE_DIR = "foam/" +TIMES = get_foam_times(CASE_DIR, string=True)[1:] + + +class TestFoamBCMappedInlet(unittest.TestCase): + """Test class for mapped inlet BCs in Hippo.""" + + def test_mapped_inlet(self): + """Test case for mapped inlet.""" + + for i in range(len(TIMES)): + self._check_u_temp_refs(i, "left", [1.0, 0, 0]) + self._check_u_temp_refs(i, "bottom", [0, 1.0, 0]) + self._check_u_temp_refs(i, "front", [0, 0, 1.0]) + + def test_mapped_inlet_subtract(self): + """Test case for mapped inlet when temperature is scaled by subtracting the difference in bulk.""" + + for i in range(len(TIMES)): + self._check_u_temp_refs(i, "left", [1.0, 0, 0], False) + self._check_u_temp_refs(i, "bottom", [0, 1.0, 0], False) + self._check_u_temp_refs(i, "front", [0, 0, 1.0], False) + + def test_mapped_inlet_rotated(self): + """Test case for when inlet's are not aligned with the axis.""" + for i in range(len(TIMES)): + self._check_u_temp_refs(i, "left", [np.sqrt(0.5), np.sqrt(0.5), 0]) + + def test_mapped_inlet_face_point(self): + """Test case for mapped inlet where the point is on the interface between boundaries.""" + + for i in range(len(TIMES)): + self._check_u_temp_refs(i, "left", [7.0 / 8, 0, 0]) + self._check_u_temp_refs(i, "bottom", [0, 7.0 / 8, 0]) + self._check_u_temp_refs(i, "front", [0, 0, 7.0 / 8]) + + def _check_u_temp_refs(self, idx, boundary, offset, use_scale=True): + rho = 0.5 + mdot_pp = 1 + t_pp = 1 + time = TIMES[idx] + + x, y, z = ff.readof.readmesh(CASE_DIR, boundary=boundary) + u = ff.readof.readvector(CASE_DIR, time, "U", boundary=boundary).T + temp = ff.readof.readscalar(CASE_DIR, time, "T", boundary=boundary) + + x += offset[0] + y += offset[1] + z += offset[2] + + if idx != 0: + t = np.float64(TIMES[idx - 1]) + u_ref = np.array( + [ + x + y + z + t, + x - y + z + t, + x + y - z + t, + ] + ).T + temp_ref = np.sqrt(x * x + y * y + z * z) + t + else: + # first time step uses initialised value + u_ref = np.array([[1, -0.5, 0.25]]) + temp_ref = 2 + + area = 4.0 + normal = np.array(offset) / np.linalg.norm(offset) + mdot = rho * np.mean(np.vecdot(u_ref, normal)) * area + u_ref *= mdot_pp / mdot + + t_bulk = np.mean(temp_ref) + if use_scale: + temp_ref *= t_pp / t_bulk + else: + temp_ref += t_pp - t_bulk + + assert np.allclose(u_ref, u, rtol=1e-7, atol=1e-12), ( + f"Max diff ({boundary}) (velocity) ({TIMES[idx]}): {abs(u - u_ref).max()} " + ) + + assert np.allclose(temp_ref, temp, rtol=1e-7, atol=1e-12), ( + f"Max diff ({boundary}) (temperature) ({time}): {abs(temp - temp_ref).max()} {temp} {temp_ref}" + ) diff --git a/test/tests/bcs/mapped_inlet/tests b/test/tests/bcs/mapped_inlet/tests new file mode 100644 index 00000000..cc19eb5a --- /dev/null +++ b/test/tests/bcs/mapped_inlet/tests @@ -0,0 +1,153 @@ +[Tests] + [mapped_mass_flow_receiver] + [setup] + type = RunCommand + command = 'bash -c "foamCleanCase -case foam && blockMesh -case foam && decomposePar -case foam"' + [] + [run_err] + type = RunException + input = main.i + prereq = mapped_mass_flow_receiver/setup + allow_warnings = true + min_parallel = 5 + max_parallel = 5 + cli_args = "FoamBCs/mass_flowx/translation_vector='-0.5 0 0'" + expect_err = "does not have a mapped plane location" + [] + [run] + type = RunApp + input = main.i + prereq = mapped_mass_flow_receiver/setup + allow_warnings = true + min_parallel = 5 + max_parallel = 5 + [] + [reconstruct] + type = RunCommand + command = 'reconstructPar -case foam' + prereq = mapped_mass_flow_receiver/run + [] + [verify] + type = PythonUnitTest + input = test.py + prereq = mapped_mass_flow_receiver/reconstruct + test_case = 'TestFoamBCMappedInlet.test_mapped_inlet' + [] + [] + [mapped_mass_flow_subtract] + [setup] + type = RunCommand + command = 'bash -c "foamCleanCase -case foam && blockMesh -case foam && decomposePar -case foam"' + [] + [run] + type = RunApp + input = main.i + prereq = mapped_mass_flow_receiver/setup + allow_warnings = true + min_parallel = 5 + max_parallel = 5 + cli_args = "FoamBCs/tempx/scale_method=subtract FoamBCs/tempy/scale_method=subtract FoamBCs/tempz/scale_method=subtract" + [] + [reconstruct] + type = RunCommand + command = 'reconstructPar -case foam' + prereq = mapped_mass_flow_receiver/run + [] + [verify] + type = PythonUnitTest + input = test.py + prereq = mapped_mass_flow_receiver/reconstruct + test_case = 'TestFoamBCMappedInlet.test_mapped_inlet_subtract' + [] + [] + [mapped_mass_flow] + [setup] + type = RunCommand + command = 'bash -c "foamCleanCase -case foam && blockMesh -case foam && decomposePar -case foam"' + prereq = mapped_mass_flow_receiver/reconstruct + [] + [run] + type = RunApp + input = main.i + prereq = mapped_mass_flow/setup + allow_warnings = true + min_parallel = 5 + max_parallel = 5 + cli_args = 'FoamBCs/mass_flowx/pp=pp FoamBCs/tempx/pp=pp FoamBCs/mass_flowy/pp=pp FoamBCs/tempy/pp=pp FoamBCs/mass_flowz/pp=pp FoamBCs/tempz/pp=pp Postprocessors/pp/expression=1' + [] + [reconstruct] + type = RunCommand + command = 'reconstructPar -case foam' + prereq = mapped_mass_flow/run + [] + [verify] + type = PythonUnitTest + input = test.py + prereq = mapped_mass_flow/reconstruct + test_case = 'TestFoamBCMappedInlet.test_mapped_inlet' + [] + [] + [mapped_rotated] + [setup] + type = RunCommand + command = 'bash -c "foamCleanCase -case foam && blockMesh -case foam && transformPoints -case foam "Rz=45" && decomposePar -case foam"' + prereq = mapped_mass_flow/reconstruct + [] + [run] + type = RunApp + input = main_rotated.i + prereq = mapped_rotated/setup + allow_warnings = true + min_parallel = 5 + max_parallel = 5 + [] + [reconstruct] + type = RunCommand + command = 'reconstructPar -case foam' + prereq = mapped_rotated/run + [] + [verify] + type = PythonUnitTest + input = test.py + prereq = mapped_rotated/reconstruct + test_case = 'TestFoamBCMappedInlet.test_mapped_inlet_rotated' + [] + [] + [face_edge_case] + [update_mesh] + type = RunCommand + command = 'bash -c "foamCleanCase -case foam && cp foam/system/blockMeshDict.face_edge foam/system/blockMeshDict && cp foam/system/decomposeParDict.face_edge foam/system/decomposeParDict"' + prereq = mapped_rotated/verify + [] + [setup] + type = RunCommand + command = 'bash -c "blockMesh -case foam && decomposePar -case foam"' + prereq = face_edge_case/update_mesh + [] + [run] + type = RunApp + input = main.i + prereq = face_edge_case/setup + allow_warnings = true + min_parallel = 2 + max_parallel = 2 + [] + [reconstruct] + type = RunCommand + command = 'reconstructPar -case foam' + prereq = face_edge_case/run + [] + [verify] + type = PythonUnitTest + input = test.py + prereq = face_edge_case/reconstruct + test_case = 'TestFoamBCMappedInlet.test_mapped_inlet_face_point' + [] + [reset_mesh] + type = RunCommand + command = 'bash -c "cp foam/system/blockMeshDict.orig foam/system/blockMeshDict && cp foam/system/decomposeParDict.orig foam/system/decomposeParDict"' + prereq = 'face_edge_case/update_mesh' + [] + [] + +[] From a3d6c20e31c00c8193e071f659a23b22e4560439 Mon Sep 17 00:00:00 2001 From: Matthew Falcone Date: Mon, 20 Apr 2026 13:40:40 +0100 Subject: [PATCH 2/7] Remove FoamControlledTimeStepper after mapped inlet rebase --- test/tests/bcs/mapped_inlet/main.i | 2 +- test/tests/bcs/mapped_inlet/main_rotated.i | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/tests/bcs/mapped_inlet/main.i b/test/tests/bcs/mapped_inlet/main.i index 50a34648..f23ae836 100644 --- a/test/tests/bcs/mapped_inlet/main.i +++ b/test/tests/bcs/mapped_inlet/main.i @@ -65,7 +65,7 @@ end_time = 5 [TimeSteppers] [foam] - type = FoamControlledTimeStepper + type = FoamTimeStepper [] [] [] diff --git a/test/tests/bcs/mapped_inlet/main_rotated.i b/test/tests/bcs/mapped_inlet/main_rotated.i index 553b6d3d..ce538978 100644 --- a/test/tests/bcs/mapped_inlet/main_rotated.i +++ b/test/tests/bcs/mapped_inlet/main_rotated.i @@ -39,7 +39,7 @@ end_time = 5 [TimeSteppers] [foam] - type = FoamControlledTimeStepper + type = FoamTimeStepper [] [] [] From c0c0448fb7b104962e02b896618fb3d0e03f4f43 Mon Sep 17 00:00:00 2001 From: Matthew Falcone Date: Mon, 8 Jun 2026 15:09:28 +0100 Subject: [PATCH 3/7] Update tests after rebase --- test/tests/bcs/mapped_inlet/tests | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/tests/bcs/mapped_inlet/tests b/test/tests/bcs/mapped_inlet/tests index cc19eb5a..234655b6 100644 --- a/test/tests/bcs/mapped_inlet/tests +++ b/test/tests/bcs/mapped_inlet/tests @@ -73,7 +73,7 @@ allow_warnings = true min_parallel = 5 max_parallel = 5 - cli_args = 'FoamBCs/mass_flowx/pp=pp FoamBCs/tempx/pp=pp FoamBCs/mass_flowy/pp=pp FoamBCs/tempy/pp=pp FoamBCs/mass_flowz/pp=pp FoamBCs/tempz/pp=pp Postprocessors/pp/expression=1' + cli_args = 'FoamBCs/mass_flowx/pp_name=pp FoamBCs/tempx/pp_name=pp FoamBCs/mass_flowy/pp_name=pp FoamBCs/tempy/pp_name=pp FoamBCs/mass_flowz/pp_name=pp FoamBCs/tempz/pp_name=pp Postprocessors/pp/expression=1' [] [reconstruct] type = RunCommand From b010c604d34f3f1c781d37bf4a48ea8fda01b82f Mon Sep 17 00:00:00 2001 From: Matthew Falcone Date: Wed, 26 Aug 2026 17:16:38 +0100 Subject: [PATCH 4/7] Update mapped inlet BC for OpenFOAM-14 --- include/base/foam/fvCFD_moose.h | 3 + include/bcs/FoamMappedInletBCBase.h | 2 + src/bcs/FoamMappedInletBCBase.C | 27 +++---- src/bcs/FoamMassFlowRateMappedInletBC.C | 3 +- src/bcs/FoamScalarBulkMappedInletBC.C | 3 +- .../mappedInletTestSolver/Make/options | 5 +- .../mappedInletTestSolver.C | 72 +------------------ .../mappedInletTestSolver.H | 43 +---------- .../bcs/mapped_inlet/foam/system/controlDict | 2 +- 9 files changed, 30 insertions(+), 130 deletions(-) diff --git a/include/base/foam/fvCFD_moose.h b/include/base/foam/fvCFD_moose.h index 61805f92..0573a4da 100644 --- a/include/base/foam/fvCFD_moose.h +++ b/include/base/foam/fvCFD_moose.h @@ -49,5 +49,8 @@ #include #include +// FoamMappedInletBCBase.h +#include + #undef NotImplemented #undef FunctionName diff --git a/include/bcs/FoamMappedInletBCBase.h b/include/bcs/FoamMappedInletBCBase.h index 29e31355..4d4cd359 100644 --- a/include/bcs/FoamMappedInletBCBase.h +++ b/include/bcs/FoamMappedInletBCBase.h @@ -23,6 +23,8 @@ class FoamMappedInletBCBase : public FoamPostprocessorBCBase MPI_Comm _mpi_comm; + const Foam::meshSearch & _mesh_searcher; + // create send and receive information for mapping void createPatchProcMap(); diff --git a/src/bcs/FoamMappedInletBCBase.C b/src/bcs/FoamMappedInletBCBase.C index 9ca25de2..7a8a67cf 100644 --- a/src/bcs/FoamMappedInletBCBase.C +++ b/src/bcs/FoamMappedInletBCBase.C @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -145,12 +146,11 @@ FoamMappedInletBCBase::createMapComm(const Foam::fvMesh & mesh, void FoamMappedInletBCBase::createPatchProcMap() { - auto & foam_mesh = _mesh->fvMesh(); - auto & boundary = foam_mesh.boundary()[_boundary[0]]; + auto & boundary = getFvMesh().boundary()[_boundary[0]]; auto face_centres = boundary.Cf(); std::vector map_procs, inlet_procs; - createMapComm(foam_mesh, face_centres, map_procs, inlet_procs); + createMapComm(getFvMesh(), face_centres, map_procs, inlet_procs); if (_mpi_comm == MPI_COMM_NULL) // process not in mapped or inlet planes return; @@ -248,16 +248,17 @@ FoamMappedInletBCBase::findIndex(const Foam::point & location, const MPI_Comm & Note that it is possible there is some non-deterministic behaviour in the function but this shouldn't be a problem in practice. */ - int index = _mesh->fvMesh().findCell(location, Foam::polyMesh::FACE_PLANES); + int index = _mesh_searcher.findCell(location, Foam::pointInCellShapes::facePlanes); int gl_index; MPI_Allreduce(&index, &gl_index, 1, MPI_INT, MPI_MAX, comm); // expand cell bounding box and repeat search if (gl_index < 0) { - Foam::label celli = _mesh->fvMesh().findNearestCell(location); + Foam::label celli = _mesh_searcher.findNearestCell(location); - bool in_cell = _mesh->fvMesh().pointInCellBB(location, celli, 0.1); + bool in_cell = + Foam::pointInCell(location, getFvMesh(), celli, Foam::pointInCellShapes::facePlanes); index = (in_cell) ? celli : -1; MPI_Allreduce(&index, &gl_index, 1, MPI_INT, MPI_MAX, comm); @@ -278,7 +279,7 @@ FoamMappedInletBCBase::findIndex(const Foam::point & location, const MPI_Comm & // use cell with cell centre closest to the location Foam::scalar dist{DBL_MAX}, gl_dist; if (index != -1) - dist = Foam::mag(_mesh->fvMesh().cellCentres()[index] - location); + dist = Foam::mag(getFvMesh().cellCentres()[index] - location); MPI_Allreduce(&dist, &gl_dist, 1, MPI_DOUBLE, MPI_MIN, comm); if (dist != gl_dist) @@ -290,7 +291,7 @@ FoamMappedInletBCBase::findIndex(const Foam::point & location, const MPI_Comm & if (in_cell > 1) { if (index != -1) - dist = Foam::mag(_mesh->fvMesh().cellCentres()[index] - (location - _offset)); + dist = Foam::mag(getFvMesh().cellCentres()[index] - (location - _offset)); MPI_Allreduce(&dist, &gl_dist, 1, MPI_DOUBLE, MPI_MIN, comm); if (dist != gl_dist) index = -1; @@ -320,12 +321,13 @@ FoamMappedInletBCBase::validParams() } FoamMappedInletBCBase::FoamMappedInletBCBase(const InputParameters & params) - : FoamPostprocessorBCBase(params), + : FoamPostprocessorBCBase(params, FoamBCType::fixedValue), _offset(), _send_map(), _recv_map(), _foam_comm(0), - _mpi_comm(MPI_COMM_NULL) + _mpi_comm(MPI_COMM_NULL), + _mesh_searcher(Foam::meshSearch::New(getFvMesh())) { if (_boundary.size() > 1) mooseError("There can only be one boundary using this method"); @@ -344,14 +346,13 @@ FoamMappedInletBCBase::getMappedArray(const Foam::word & name) if (_mpi_comm == MPI_COMM_NULL) return Foam::Field(); - auto & foam_mesh = _mesh->fvMesh(); - auto & boundary_patch = foam_mesh.boundary()[_boundary[0]]; + auto & boundary_patch = getFvMesh().boundary()[_boundary[0]]; Foam::PstreamBuffers sendBuf( Foam::UPstream::commsTypes::nonBlocking, Foam::UPstream::msgType(), _foam_comm); if (_send_map.size() > 0) { - auto & var = foam_mesh.lookupObject>(name); + auto & var = getFvMesh().lookupObject>(name); for (auto & pair : _send_map) { auto & proc = pair.first; diff --git a/src/bcs/FoamMassFlowRateMappedInletBC.C b/src/bcs/FoamMassFlowRateMappedInletBC.C index 88ba1f96..60a0e499 100644 --- a/src/bcs/FoamMassFlowRateMappedInletBC.C +++ b/src/bcs/FoamMassFlowRateMappedInletBC.C @@ -39,8 +39,7 @@ FoamMassFlowRateMappedInletBC::FoamMassFlowRateMappedInletBC(const InputParamete void FoamMassFlowRateMappedInletBC::imposeBoundaryCondition() { - auto & foam_mesh = _mesh->fvMesh(); - auto & boundary_patch = foam_mesh.boundary()[_boundary[0]]; + auto & boundary_patch = getFvMesh().boundary()[_boundary[0]]; // currently we map mass flux rather than velocity, maybe useful to have option auto && U_map = getMappedArray("U"); diff --git a/src/bcs/FoamScalarBulkMappedInletBC.C b/src/bcs/FoamScalarBulkMappedInletBC.C index f1cb83d9..c0d34701 100644 --- a/src/bcs/FoamScalarBulkMappedInletBC.C +++ b/src/bcs/FoamScalarBulkMappedInletBC.C @@ -32,8 +32,7 @@ FoamScalarBulkMappedInletBC::FoamScalarBulkMappedInletBC(const InputParameters & void FoamScalarBulkMappedInletBC::imposeBoundaryCondition() { - auto & foam_mesh = _mesh->fvMesh(); - auto & boundary_patch = foam_mesh.boundary()[_boundary[0]]; + auto & boundary_patch = getFvMesh().boundary()[_boundary[0]]; auto && var_map = getMappedArray(_foam_variable); auto & Sf = boundary_patch.magSf(); diff --git a/test/OpenFOAM/modules/mappedInletTestSolver/Make/options b/test/OpenFOAM/modules/mappedInletTestSolver/Make/options index cc845ca8..37d3496e 100644 --- a/test/OpenFOAM/modules/mappedInletTestSolver/Make/options +++ b/test/OpenFOAM/modules/mappedInletTestSolver/Make/options @@ -6,9 +6,12 @@ EXE_INC = \ -I$(LIB_SRC)/ThermophysicalTransportModels/solid/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(LIB_SRC)/sampling/lnInclude + -I$(LIB_SRC)/sampling/lnInclude \ + -I../baseTestSolver LIB_LIBS = \ + -L$(FOAM_USER_LIBBIN) \ + -lbaseTestSolver \ -lsolidThermo \ -lsolidThermophysicalTransportModels \ -lcoupledThermophysicalTransportModels \ diff --git a/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.C b/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.C index 17f910c8..1afdfec6 100644 --- a/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.C +++ b/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.C @@ -23,7 +23,6 @@ License \*---------------------------------------------------------------------------*/ -#include "dimensionSets.H" #include "dimensionedScalar.H" #include "dimensionedVector.H" #include "fvMesh.H" @@ -48,29 +47,11 @@ addToRunTimeSelectionTable(solver, mappedInletTestSolver, fvMesh); // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // -bool -Foam::solvers::mappedInletTestSolver::dependenciesModified() const -{ - return runTime.controlDict().modified(); -} - -bool -Foam::solvers::mappedInletTestSolver::read() -{ - solver::read(); - - maxDeltaT_ = runTime.controlDict().found("maxDeltaT") - ? runTime.controlDict().lookup("maxDeltaT", runTime.userUnits()) - : vGreat; - - return true; -} - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // Solver based on solid.C module Foam::solvers::mappedInletTestSolver::mappedInletTestSolver(fvMesh & mesh, autoPtr thermoPtr) - : solver(mesh), + : baseTestSolver(mesh), thermoPtr_(thermoPtr), thermo_(thermoPtr_()), @@ -94,18 +75,8 @@ Foam::solvers::mappedInletTestSolver::mappedInletTestSolver(fvMesh & mesh) read(); } -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::solvers::mappedInletTestSolver::~mappedInletTestSolver() {} - // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // -Foam::scalar -Foam::solvers::mappedInletTestSolver::maxDeltaT() const -{ - return min(fvModels().maxDeltaT(), maxDeltaT_); -} - void Foam::solvers::mappedInletTestSolver::preSolve() { @@ -114,32 +85,6 @@ Foam::solvers::mappedInletTestSolver::preSolve() // Update the mesh for topology change, mesh to mesh mapping mesh_.update(); } - -void -Foam::solvers::mappedInletTestSolver::moveMesh() -{ - if (pimple.firstIter() || pimple.moveMeshOuterCorrectors()) - { - if (!mesh_.mover().solidBody()) - { - FatalErrorInFunction << "Region " << name() << " of type " << type() - << " does not support non-solid body mesh motion" << exit(FatalError); - } - - mesh_.move(); - } -} - -void -Foam::solvers::mappedInletTestSolver::motionCorrector() -{ -} - -void -Foam::solvers::mappedInletTestSolver::prePredictor() -{ -} - void Foam::solvers::mappedInletTestSolver::momentumPredictor() { @@ -165,19 +110,4 @@ Foam::solvers::mappedInletTestSolver::thermophysicalPredictor() thermo_.correct(); } -void -Foam::solvers::mappedInletTestSolver::pressureCorrector() -{ -} - -void -Foam::solvers::mappedInletTestSolver::postCorrector() -{ -} - -void -Foam::solvers::mappedInletTestSolver::postSolve() -{ -} - // ************************************************************************* // diff --git a/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.H b/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.H index 24705868..a64d4ce7 100644 --- a/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.H +++ b/test/OpenFOAM/modules/mappedInletTestSolver/mappedInletTestSolver.H @@ -34,7 +34,7 @@ SourceFiles #pragma once -#include "solver.H" +#include "baseTestSolver.H" #include "solidThermophysicalTransportModel.H" #include "volFieldsFwd.H" @@ -49,7 +49,7 @@ namespace solvers Class mappedInletTestSolver Declaration \*---------------------------------------------------------------------------*/ -class mappedInletTestSolver : public solver +class mappedInletTestSolver : public baseTestSolver { protected: @@ -67,14 +67,6 @@ protected: autoPtr thermophysicalTransport; - // Protected Member Functions - - //- Return true if the solver's dependencies have been modified - virtual bool dependenciesModified() const; - - //- Read controls - virtual bool read(); - public: // Public Data // reference to thermophysical properties @@ -87,7 +79,7 @@ public: const volVectorField & U; //- Runtime type information - TypeName("mappedInletTestSolver"); + TypeName("mappedInletTest"); // Constructors @@ -99,26 +91,11 @@ public: //- Disallow default bitwise copy construction mappedInletTestSolver(const mappedInletTestSolver &) = delete; - //- Destructor - virtual ~mappedInletTestSolver(); - // Member Functions - //- Return the current maximum time-step for stable solution - virtual scalar maxDeltaT() const; - //- Called at the start of the time-step, before the PIMPLE loop virtual void preSolve(); - //- Called at the start of the PIMPLE loop to move the mesh - virtual void moveMesh(); - - //- Corrections that follow mesh motion - virtual void motionCorrector(); - - //- Called at the beginning of the PIMPLE loop - virtual void prePredictor(); - //- Construct and optionally solve the momentum equation virtual void momentumPredictor(); @@ -126,20 +103,6 @@ public: // convert to temperature // and update thermophysical and transport properties virtual void thermophysicalPredictor(); - - //- Construct and solve the pressure equation in the PISO loop - virtual void pressureCorrector(); - - //- Correct the thermophysical transport modelling - virtual void postCorrector(); - - //- Called after the PIMPLE loop at the end of the time-step - virtual void postSolve(); - - // Member Operators - - //- Disallow default bitwise assignment - void operator=(const mappedInletTestSolver &) = delete; }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/test/tests/bcs/mapped_inlet/foam/system/controlDict b/test/tests/bcs/mapped_inlet/foam/system/controlDict index 2bbe0e5f..29bdb54c 100644 --- a/test/tests/bcs/mapped_inlet/foam/system/controlDict +++ b/test/tests/bcs/mapped_inlet/foam/system/controlDict @@ -13,7 +13,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -solver mappedInletTestSolver; +solver mappedInletTest; startFrom startTime; From 338fb4cbd26ea80f58c16694a4f576569d6b9b72 Mon Sep 17 00:00:00 2001 From: Matthew Falcone Date: Thu, 27 Aug 2026 15:43:21 +0100 Subject: [PATCH 5/7] Initial pass at improving mapped inlet implementation --- include/bcs/FoamMappedInletBCBase.h | 13 +- src/bcs/FoamMappedInletBCBase.C | 311 +++++++++++++++----------- src/bcs/FoamScalarBulkMappedInletBC.C | 5 + test/tests/bcs/mapped_inlet/tests | 14 +- 4 files changed, 211 insertions(+), 132 deletions(-) diff --git a/include/bcs/FoamMappedInletBCBase.h b/include/bcs/FoamMappedInletBCBase.h index 4d4cd359..8b205119 100644 --- a/include/bcs/FoamMappedInletBCBase.h +++ b/include/bcs/FoamMappedInletBCBase.h @@ -2,6 +2,8 @@ #include "FoamPostprocessorBCBase.h" #include +#include +#include class FoamMappedInletBCBase : public FoamPostprocessorBCBase { @@ -37,12 +39,13 @@ class FoamMappedInletBCBase : public FoamPostprocessorBCBase // create/assign communicators for the transfers between map and inlet planes void createMapComm(const Foam::fvMesh & mesh, - Foam::vectorField face_centres, + const Foam::vectorField & face_centres, std::vector & send_process, - std::vector & recv_process); + std::vector & recv_process, + const MPI_Comm & comm); - // find index of cell containing point or raise error if not found - int findIndex(const Foam::point & location, const MPI_Comm & comm); + // find indices of cells containing mapped inlet points or raise error if not found + std::vector findIndices(const Foam::pointField & locations, const MPI_Comm & comm); // handle creation of new communicators in parallel or serial Foam::label @@ -58,7 +61,7 @@ class FoamMappedInletBCBase : public FoamPostprocessorBCBase else { foam_comm = Foam::UPstream::worldComm; - new_comm = MPI_COMM_WORLD; + new_comm = Foam::PstreamGlobals::MPICommunicators_[Foam::UPstream::worldComm]; } return foam_comm; } diff --git a/src/bcs/FoamMappedInletBCBase.C b/src/bcs/FoamMappedInletBCBase.C index 7a8a67cf..c9c170db 100644 --- a/src/bcs/FoamMappedInletBCBase.C +++ b/src/bcs/FoamMappedInletBCBase.C @@ -1,30 +1,34 @@ #include "FoamMappedInletBCBase.h" +#include "MooseError.h" #include #include #include #include #include +#include +#include +#include #include #include #include #include #include -#include +#include namespace { // Get the cartesian bounding box of the mapped inlet plane void -getBBox(const Foam::vectorField points, Real bbox[6]) +getBBox(const Foam::vectorField & points, Real bbox[6], const MPI_Comm & comm) { - bbox[0] = DBL_MAX; - bbox[1] = DBL_MIN; - bbox[2] = DBL_MAX; - bbox[3] = DBL_MIN; - bbox[4] = DBL_MAX; - bbox[5] = DBL_MIN; + bbox[0] = std::numeric_limits::max(); + bbox[1] = std::numeric_limits::lowest(); + bbox[2] = std::numeric_limits::max(); + bbox[3] = std::numeric_limits::lowest(); + bbox[4] = std::numeric_limits::max(); + bbox[5] = std::numeric_limits::lowest(); for (auto p : points) { bbox[0] = std::min(bbox[0], p.x()); @@ -35,12 +39,42 @@ getBBox(const Foam::vectorField points, Real bbox[6]) bbox[5] = std::max(bbox[5], p.z()); } - MPI_Allreduce(MPI_IN_PLACE, &bbox[0], 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, &bbox[1], 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, &bbox[2], 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, &bbox[3], 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, &bbox[4], 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, &bbox[5], 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &bbox[0], 1, MPI_DOUBLE, MPI_MIN, comm); + MPI_Allreduce(MPI_IN_PLACE, &bbox[1], 1, MPI_DOUBLE, MPI_MAX, comm); + MPI_Allreduce(MPI_IN_PLACE, &bbox[2], 1, MPI_DOUBLE, MPI_MIN, comm); + MPI_Allreduce(MPI_IN_PLACE, &bbox[3], 1, MPI_DOUBLE, MPI_MAX, comm); + MPI_Allreduce(MPI_IN_PLACE, &bbox[4], 1, MPI_DOUBLE, MPI_MIN, comm); + MPI_Allreduce(MPI_IN_PLACE, &bbox[5], 1, MPI_DOUBLE, MPI_MAX, comm); +} + +enum class FindError +{ + NOT_FOUND, + MANY_FOUND +}; + +std::vector +findLocationErrors(const std::vector & indices, const MPI_Comm & comm, FindError issue) +{ + std::vector found(indices.size(), 0); + std::vector gl_found(indices.size(), 0); + + for (auto i = 0lu; i < indices.size(); ++i) + { + found[i] = indices[i] >= 0; + } + MPI_Allreduce(found.data(), gl_found.data(), found.size(), MPI_INT, MPI_SUM, comm); + + std::vector deconflict_points; + for (auto i = 0lu; i < indices.size(); ++i) + { + if (issue == FindError::MANY_FOUND && gl_found[i] > 1) + deconflict_points.push_back(i); + else if (issue == FindError::NOT_FOUND && gl_found[i] == 0) + deconflict_points.push_back((i)); + } + + return deconflict_points; } } @@ -51,20 +85,12 @@ FoamMappedInletBCBase::intersectMapPlane(const Foam::fvMesh & mesh, Real cart_bb for (int i = 0; i < mesh.nCells(); ++i) { auto points = mesh.cellPoints(i); - for (auto pointI : points) - { - // check whether cell intersects bbox - const Foam::point & p = vertices[pointI]; - if (p.x() >= cart_bbox[0] && p.x() <= cart_bbox[1] && p.y() >= cart_bbox[2] && - p.y() <= cart_bbox[3] && p.z() >= cart_bbox[4] && p.z() <= cart_bbox[5]) - { - return true; - } - } - - // The bbox could by narrower than the cell, check in each direction whether bbox is within the - // cell bbox - Real cell_bbox[6] = {DBL_MAX, DBL_MIN, DBL_MAX, DBL_MIN, DBL_MAX, DBL_MIN}; + Real cell_bbox[6] = {std::numeric_limits::max(), + std::numeric_limits::lowest(), + std::numeric_limits::max(), + std::numeric_limits::lowest(), + std::numeric_limits::max(), + std::numeric_limits::lowest()}; for (auto point : points) { const Foam::point & p = vertices[point]; @@ -75,25 +101,9 @@ FoamMappedInletBCBase::intersectMapPlane(const Foam::fvMesh & mesh, Real cart_bb cell_bbox[4] = std::min(cell_bbox[4], p.z()); cell_bbox[5] = std::max(cell_bbox[5], p.z()); } - // check is cart bbox is narrower than cell bbox in x direction - if ((cart_bbox[0] >= cell_bbox[0] && cart_bbox[1] <= cell_bbox[1]) && - (cart_bbox[3] >= cell_bbox[2] || cart_bbox[2] <= cell_bbox[3]) && - (cart_bbox[5] >= cell_bbox[4] || cart_bbox[4] <= cell_bbox[5])) - { - return true; - } - - // check is cart bbox is narrower than cell bbox in y direction - if ((cart_bbox[2] >= cell_bbox[2] && cart_bbox[3] <= cell_bbox[3]) && - (cart_bbox[1] >= cell_bbox[0] || cart_bbox[0] <= cell_bbox[1]) && - (cart_bbox[5] >= cell_bbox[4] || cart_bbox[4] <= cell_bbox[5])) - { - return true; - } - // check is cart bbox is narrower than cell bbox in z direction - if ((cart_bbox[4] >= cell_bbox[4] && cart_bbox[5] <= cell_bbox[5]) && - (cart_bbox[1] >= cell_bbox[0] || cart_bbox[0] <= cell_bbox[1]) && - (cart_bbox[3] >= cell_bbox[2] || cart_bbox[2] <= cell_bbox[3])) + if (cart_bbox[0] <= cell_bbox[1] && cart_bbox[1] >= cell_bbox[0] && + cart_bbox[2] <= cell_bbox[3] && cart_bbox[3] >= cell_bbox[2] && + cart_bbox[4] <= cell_bbox[5] && cart_bbox[5] >= cell_bbox[4]) { return true; } @@ -104,13 +114,14 @@ FoamMappedInletBCBase::intersectMapPlane(const Foam::fvMesh & mesh, Real cart_bb void FoamMappedInletBCBase::createMapComm(const Foam::fvMesh & mesh, - const Foam::vectorField face_centres, + const Foam::vectorField & face_centres, std::vector & map_process, - std::vector & inlet_process) + std::vector & inlet_process, + const MPI_Comm & comm) { Real cart_bbox[6]; auto mapped_plane = face_centres + _offset; - getBBox(mapped_plane(), cart_bbox); + getBBox(mapped_plane(), cart_bbox, comm); int mappedPlaneProcess = intersectMapPlane(mesh, cart_bbox); int inletPlaneProcess = face_centres.size() > 0; @@ -118,8 +129,8 @@ FoamMappedInletBCBase::createMapComm(const Foam::fvMesh & mesh, std::vector inlet_procs(Foam::UPstream::nProcs()); std::vector map_procs(Foam::UPstream::nProcs()); - MPI_Allgather(&mappedPlaneProcess, 1, MPI_INT, map_procs.data(), 1, MPI_INT, MPI_COMM_WORLD); - MPI_Allgather(&inletPlaneProcess, 1, MPI_INT, inlet_procs.data(), 1, MPI_INT, MPI_COMM_WORLD); + MPI_Allgather(&mappedPlaneProcess, 1, MPI_INT, map_procs.data(), 1, MPI_INT, comm); + MPI_Allgather(&inletPlaneProcess, 1, MPI_INT, inlet_procs.data(), 1, MPI_INT, comm); map_process.clear(); inlet_process.clear(); @@ -150,7 +161,11 @@ FoamMappedInletBCBase::createPatchProcMap() auto face_centres = boundary.Cf(); std::vector map_procs, inlet_procs; - createMapComm(getFvMesh(), face_centres, map_procs, inlet_procs); + const MPI_Comm & comm{Foam::PstreamGlobals::MPICommunicators_[Foam::UPstream::worldComm]}; + createMapComm(getFvMesh(), face_centres, map_procs, inlet_procs, comm); + + if (map_procs.empty()) + mooseError("Mapped inlet invalid. No intersection of mapped plane and mesh."); if (_mpi_comm == MPI_COMM_NULL) // process not in mapped or inlet planes return; @@ -172,6 +187,7 @@ FoamMappedInletBCBase::createPatchProcMap() send_points.finishedSends(true); + std::vector send_sizes(inlet_procs.size()); std::vector size_requests(inlet_procs.size()); std::vector data_requests(inlet_procs.size()); std::vector> recv_indices_procs(inlet_procs.size()); @@ -187,12 +203,13 @@ FoamMappedInletBCBase::createPatchProcMap() recieve >> field; auto & vec = _send_map[inlet_procs[i]]; auto & recv_indices = recv_indices_procs[i]; - for (int j = 0; j < field.size(); ++j) + auto indices = findIndices(field + _offset, map_comm); + + for (auto j = 0lu; j < indices.size(); ++j) { - auto index = findIndex(field[j] + _offset, map_comm); - if (index >= 0) + if (indices[j] >= 0) { - vec.push_back(index); // assign to send map required indices + vec.push_back(indices[j]); // assign to send map required indices recv_indices.push_back(j); } } @@ -200,8 +217,8 @@ FoamMappedInletBCBase::createPatchProcMap() _send_map.erase(inlet_procs[i]); // Let original processes know which points will come from each rank - int size = recv_indices.size(); - MPI_Isend(&size, 1, MPI_INT, inlet_procs[i], 0, _mpi_comm, &size_requests.at(i)); + send_sizes[i] = recv_indices.size(); + MPI_Isend(&send_sizes[i], 1, MPI_INT, inlet_procs[i], 0, _mpi_comm, &size_requests.at(i)); MPI_Isend(recv_indices.data(), recv_indices.size(), MPI_INT, @@ -216,98 +233,142 @@ FoamMappedInletBCBase::createPatchProcMap() if (isInletProc) // create map to determine where data from map processes should go { + std::vector assignment_count(face_centres.size(), 0); for (auto & proc : map_procs) { - int size; + int size{0}; MPI_Recv(&size, 1, MPI_INT, proc, 0, _mpi_comm, MPI_STATUS_IGNORE); std::vector recv_indices(size); MPI_Recv(recv_indices.data(), size, MPI_INT, proc, 1, _mpi_comm, MPI_STATUS_IGNORE); for (auto & index : recv_indices) { - assert(index < face_centres.size()); + if (index < 0 || index >= face_centres.size()) + mooseError( + "Received invalid mapped inlet face index ", index, " from process ", proc, "."); + _recv_map[proc].push_back(index); + ++assignment_count[index]; } } + + for (auto i = 0lu; i < assignment_count.size(); ++i) + if (assignment_count[i] != 1) + mooseError("Mapped inlet face at (", + face_centres[i].x(), + ", ", + face_centres[i].y(), + ", ", + face_centres[i].z(), + ") was assigned to ", + assignment_count[i], + " mapped cells; expected exactly one."); + } + + if (isMapProc) + { + MPI_Waitall(inlet_procs.size(), size_requests.data(), MPI_STATUSES_IGNORE); + MPI_Waitall(inlet_procs.size(), data_requests.data(), MPI_STATUSES_IGNORE); } - MPI_Barrier(_mpi_comm); } -int -FoamMappedInletBCBase::findIndex(const Foam::point & location, const MPI_Comm & comm) +std::vector +FoamMappedInletBCBase::findIndices(const Foam::pointField & locations, const MPI_Comm & comm) { - /* - This function uses several ways of finding the mapped plane cell - 1. use findCell function - - Sometimes on cell boundaries it may not find the cell or two processes with both find it - 2. If none found, find the closest point and do an expanded bounding box search, raise error - if none still found. - 3. If multiple found, use the cell with cell centre closer to desired location - 4. If still multiple found, use the cell closer to the inlet. - - In the unlikely chance there are still multiple cells detected, raise a warning - Note that it is possible there is some non-deterministic behaviour in the function but this - shouldn't be a problem in practice. - */ - int index = _mesh_searcher.findCell(location, Foam::pointInCellShapes::facePlanes); - int gl_index; - MPI_Allreduce(&index, &gl_index, 1, MPI_INT, MPI_MAX, comm); - - // expand cell bounding box and repeat search - if (gl_index < 0) + std::vector indices(locations.size(), -1); + for (auto i = 0; i < locations.size(); ++i) { - Foam::label celli = _mesh_searcher.findNearestCell(location); + indices[i] = _mesh_searcher.findCell(locations[i], Foam::pointInCellShapes::facePlanes); + } - bool in_cell = - Foam::pointInCell(location, getFvMesh(), celli, Foam::pointInCellShapes::facePlanes); - index = (in_cell) ? celli : -1; - MPI_Allreduce(&index, &gl_index, 1, MPI_INT, MPI_MAX, comm); + auto bad_indices = findLocationErrors(indices, comm, FindError::NOT_FOUND); + for (auto i = 0lu; i < bad_indices.size(); ++i) + { + int celli = _mesh_searcher.findNearestCell(locations[bad_indices[i]]); + bool in_cell = Foam::pointInCell( + locations[bad_indices[i]], getFvMesh(), celli, Foam::pointInCellShapes::facePlanes); + indices[bad_indices[i]] = (in_cell) ? celli : -1; + } - int rank; - MPI_Comm_rank(comm, &rank); - if (gl_index < 0 && rank == 0) + // Check not founds + std::vector notfound_points = findLocationErrors(indices, comm, FindError::NOT_FOUND); + if (notfound_points.size() > 0) + { + std::stringstream err_msg; + for (auto idx : notfound_points) { - mooseError("Face centre at location (", - location[0], - ",", - location[1], - ",", - location[2], - ") does not have a mapped plane location"); + const auto & loc = locations[idx]; + err_msg << "\t(" << loc.x() << ", " << loc.y() << ", " << loc.z() << ")\n"; } + mooseError("Locations not found:\n", err_msg.str()); } - // use cell with cell centre closest to the location - Foam::scalar dist{DBL_MAX}, gl_dist; - if (index != -1) - dist = Foam::mag(getFvMesh().cellCentres()[index] - location); + // Deduplication + // check if deduplication is required + auto deconflict_points = findLocationErrors(indices, comm, FindError::MANY_FOUND); + if (deconflict_points.empty()) + return indices; - MPI_Allreduce(&dist, &gl_dist, 1, MPI_DOUBLE, MPI_MIN, comm); - if (dist != gl_dist) - index = -1; + // deconflict step 1 + // Use cell centre closest to location + std::vector dist(deconflict_points.size(), std::numeric_limits::max()); + std::vector gl_dist(deconflict_points.size(), std::numeric_limits::max()); - // 2. use cell centre closest to inlet point - int in_cell = index != -1; - MPI_Allreduce(MPI_IN_PLACE, &in_cell, 1, MPI_INT, MPI_SUM, comm); - if (in_cell > 1) + for (auto i = 0lu; i < deconflict_points.size(); ++i) { - if (index != -1) - dist = Foam::mag(getFvMesh().cellCentres()[index] - (location - _offset)); - MPI_Allreduce(&dist, &gl_dist, 1, MPI_DOUBLE, MPI_MIN, comm); - if (dist != gl_dist) - index = -1; + int idx = deconflict_points[i]; + + if (indices[idx] != -1) + dist[i] = Foam::mag(locations[idx] - getFvMesh().cellCentres()[indices[idx]]); } - in_cell = index != -1; - MPI_Allreduce(MPI_IN_PLACE, &in_cell, 1, MPI_INT, MPI_SUM, comm); - - if (in_cell > 1) - mooseWarning("More than 1 process found location (", - location[0], - ",", - location[1], - ",", - location[2], - ")"); - return index; + + MPI_Allreduce(dist.data(), gl_dist.data(), deconflict_points.size(), MPI_DOUBLE, MPI_MIN, comm); + + for (auto i = 0lu; i < deconflict_points.size(); ++i) + { + int idx = deconflict_points[i]; + if (dist[i] != gl_dist[i]) + indices[idx] = -1; + } + + // Check again for conflicts + deconflict_points = findLocationErrors(indices, comm, FindError::MANY_FOUND); + if (deconflict_points.empty()) + return indices; + + // 2. Choose which ever is closest to the inlet + dist.assign(deconflict_points.size(), std::numeric_limits::max()); + gl_dist.assign(deconflict_points.size(), std::numeric_limits::max()); + for (auto i = 0lu; i < deconflict_points.size(); ++i) + { + int idx = deconflict_points[i]; + + if (indices[idx] != -1) + dist[i] = Foam::mag(getFvMesh().cellCentres()[indices[idx]] - (locations[idx] - _offset)); + } + + MPI_Allreduce(dist.data(), gl_dist.data(), deconflict_points.size(), MPI_DOUBLE, MPI_MIN, comm); + + for (auto i = 0lu; i < deconflict_points.size(); ++i) + { + int idx = deconflict_points[i]; + if (dist[i] != gl_dist[i]) + indices[idx] = -1; + } + + deconflict_points = findLocationErrors(indices, comm, FindError::MANY_FOUND); + if (!deconflict_points.empty()) + { + std::stringstream err_msg; + for (auto idx : deconflict_points) + { + const auto & loc = locations[idx]; + err_msg << "\t(" << loc.x() << ", " << loc.y() << ", " << loc.z() << ")\n"; + } + mooseError("Locations found in multiple ranks\n", err_msg.str()); + } + + return indices; } InputParameters diff --git a/src/bcs/FoamScalarBulkMappedInletBC.C b/src/bcs/FoamScalarBulkMappedInletBC.C index c0d34701..54ca4388 100644 --- a/src/bcs/FoamScalarBulkMappedInletBC.C +++ b/src/bcs/FoamScalarBulkMappedInletBC.C @@ -1,11 +1,13 @@ #include "FoamScalarBulkMappedInletBC.h" #include "InputParameters.h" +#include "MooseError.h" #include "MooseTypes.h" #include "PstreamReduceOps.H" #include "Registry.h" #include "ops.H" #include "volFieldsFwd.H" +#include registerMooseObject("hippoApp", FoamScalarBulkMappedInletBC); @@ -56,6 +58,9 @@ FoamScalarBulkMappedInletBC::applyScaleMethod(T & var, const Real bulk_ref, cons { if (_scale_method == "SCALE") { + if (std::abs(bulk) < Foam::SMALL) + mooseError("Cannot scale mapped inlet field '", _foam_variable, "', bulk value is 0."); + return (var * bulk_ref / bulk)(); } else if (_scale_method == "SUBTRACT") diff --git a/test/tests/bcs/mapped_inlet/tests b/test/tests/bcs/mapped_inlet/tests index 234655b6..4ef88cfa 100644 --- a/test/tests/bcs/mapped_inlet/tests +++ b/test/tests/bcs/mapped_inlet/tests @@ -4,7 +4,7 @@ type = RunCommand command = 'bash -c "foamCleanCase -case foam && blockMesh -case foam && decomposePar -case foam"' [] - [run_err] + [run_err1] type = RunException input = main.i prereq = mapped_mass_flow_receiver/setup @@ -12,7 +12,17 @@ min_parallel = 5 max_parallel = 5 cli_args = "FoamBCs/mass_flowx/translation_vector='-0.5 0 0'" - expect_err = "does not have a mapped plane location" + expect_err = "Mapped inlet invalid. No intersection of mapped plane and mesh." + [] + [run_err2] + type = RunException + input = main.i + prereq = mapped_mass_flow_receiver/setup + allow_warnings = true + min_parallel = 5 + max_parallel = 5 + cli_args = "FoamBCs/mass_flowx/translation_vector='1 0.5 0'" + expect_err = "Locations not found:" [] [run] type = RunApp From e474a5f1ac81919c36c2cafadffab2a376149b5d Mon Sep 17 00:00:00 2001 From: Matthew Falcone Date: Thu, 27 Aug 2026 16:30:16 +0100 Subject: [PATCH 6/7] Clean up and refactor findIndices --- include/bcs/FoamMappedInletBCBase.h | 37 ++---- src/bcs/FoamMappedInletBCBase.C | 177 ++++++++++++++------------ src/bcs/FoamScalarBulkMappedInletBC.C | 6 +- test/tests/bcs/mapped_inlet/tests | 2 +- test/tests/bcs/receiver_pp/test.py | 2 +- test/tests/bcs/receiver_pp/tests | 2 +- 6 files changed, 115 insertions(+), 111 deletions(-) diff --git a/include/bcs/FoamMappedInletBCBase.h b/include/bcs/FoamMappedInletBCBase.h index 8b205119..09542b1e 100644 --- a/include/bcs/FoamMappedInletBCBase.h +++ b/include/bcs/FoamMappedInletBCBase.h @@ -27,13 +27,14 @@ class FoamMappedInletBCBase : public FoamPostprocessorBCBase const Foam::meshSearch & _mesh_searcher; - // create send and receive information for mapping - void createPatchProcMap(); - // get array from mapped plane on the inlet processes template Foam::Field getMappedArray(const Foam::word & name); +private: + // create send and receive information for mapping + void createPatchProcMap(); + // check if bounding box intersects with rank bool intersectMapPlane(const Foam::fvMesh & mesh, Real cart_bbox[6]); @@ -42,34 +43,16 @@ class FoamMappedInletBCBase : public FoamPostprocessorBCBase const Foam::vectorField & face_centres, std::vector & send_process, std::vector & recv_process, - const MPI_Comm & comm); + MPI_Comm comm); // find indices of cells containing mapped inlet points or raise error if not found - std::vector findIndices(const Foam::pointField & locations, const MPI_Comm & comm); + std::vector findIndices(const Foam::pointField & locations, MPI_Comm comm); // handle creation of new communicators in parallel or serial - Foam::label - createCommunicator(const Foam::label parent_comm, std::vector procs, MPI_Comm & new_comm) - { - Foam::label foam_comm; - if (Foam::UPstream::parRun()) - { - Foam::labelList foam_procs(procs.begin(), procs.end()); - foam_comm = Foam::UPstream::allocateCommunicator(parent_comm, foam_procs, true); - new_comm = Foam::PstreamGlobals::MPICommunicators_[foam_comm]; - } - else - { - foam_comm = Foam::UPstream::worldComm; - new_comm = Foam::PstreamGlobals::MPICommunicators_[Foam::UPstream::worldComm]; - } - return foam_comm; - } + Foam::label createCommunicator(const Foam::label parent_comm, + const std::vector & procs, + MPI_Comm & new_comm); // free communicators if parallel run - void destroyCommunicator(Foam::label comm) - { - if (Foam::UPstream::parRun()) - Foam::UPstream::freeCommunicator(comm); - } + void destroyCommunicator(Foam::label comm); }; diff --git a/src/bcs/FoamMappedInletBCBase.C b/src/bcs/FoamMappedInletBCBase.C index c9c170db..2c271897 100644 --- a/src/bcs/FoamMappedInletBCBase.C +++ b/src/bcs/FoamMappedInletBCBase.C @@ -1,18 +1,17 @@ #include "FoamMappedInletBCBase.h" -#include "MooseError.h" #include #include + #include #include #include -#include -#include -#include #include #include #include +#include +#include #include #include #include @@ -21,7 +20,7 @@ namespace { // Get the cartesian bounding box of the mapped inlet plane void -getBBox(const Foam::vectorField & points, Real bbox[6], const MPI_Comm & comm) +getBBox(const Foam::vectorField & points, Real bbox[6], MPI_Comm comm) { bbox[0] = std::numeric_limits::max(); bbox[1] = std::numeric_limits::lowest(); @@ -47,14 +46,14 @@ getBBox(const Foam::vectorField & points, Real bbox[6], const MPI_Comm & comm) MPI_Allreduce(MPI_IN_PLACE, &bbox[5], 1, MPI_DOUBLE, MPI_MAX, comm); } -enum class FindError +enum class LocatorIssue { NOT_FOUND, MANY_FOUND }; std::vector -findLocationErrors(const std::vector & indices, const MPI_Comm & comm, FindError issue) +findLocatorIssues(const std::vector & indices, MPI_Comm comm, LocatorIssue issue) { std::vector found(indices.size(), 0); std::vector gl_found(indices.size(), 0); @@ -65,16 +64,45 @@ findLocationErrors(const std::vector & indices, const MPI_Comm & comm, Find } MPI_Allreduce(found.data(), gl_found.data(), found.size(), MPI_INT, MPI_SUM, comm); - std::vector deconflict_points; + std::vector conflict_indices; for (auto i = 0lu; i < indices.size(); ++i) { - if (issue == FindError::MANY_FOUND && gl_found[i] > 1) - deconflict_points.push_back(i); - else if (issue == FindError::NOT_FOUND && gl_found[i] == 0) - deconflict_points.push_back((i)); + if (issue == LocatorIssue::MANY_FOUND && gl_found[i] > 1) + conflict_indices.push_back(i); + else if (issue == LocatorIssue::NOT_FOUND && gl_found[i] == 0) + conflict_indices.push_back(i); + } + + return conflict_indices; +} + +void +retainClosestCells(std::vector & indices, + const std::vector & conflict_indices, + const Foam::pointField & target_locations, + const Foam::vectorField & cell_centres, + MPI_Comm comm) +{ + std::vector local_distances(conflict_indices.size(), std::numeric_limits::max()); + std::vector minimum_distances(conflict_indices.size()); + + for (auto i = 0lu; i < conflict_indices.size(); ++i) + { + const int idx = conflict_indices[i]; + if (indices[idx] >= 0) + local_distances[i] = Foam::mag(target_locations[idx] - cell_centres[indices[idx]]); } - return deconflict_points; + MPI_Allreduce(local_distances.data(), + minimum_distances.data(), + conflict_indices.size(), + MPI_DOUBLE, + MPI_MIN, + comm); + + for (auto i = 0lu; i < conflict_indices.size(); ++i) + if (local_distances[i] != minimum_distances[i]) + indices[conflict_indices[i]] = -1; } } @@ -117,14 +145,14 @@ FoamMappedInletBCBase::createMapComm(const Foam::fvMesh & mesh, const Foam::vectorField & face_centres, std::vector & map_process, std::vector & inlet_process, - const MPI_Comm & comm) + MPI_Comm comm) { Real cart_bbox[6]; auto mapped_plane = face_centres + _offset; getBBox(mapped_plane(), cart_bbox, comm); int mappedPlaneProcess = intersectMapPlane(mesh, cart_bbox); - int inletPlaneProcess = face_centres.size() > 0; + int inletPlaneProcess = !face_centres.empty(); std::vector inlet_procs(Foam::UPstream::nProcs()); std::vector map_procs(Foam::UPstream::nProcs()); @@ -165,7 +193,7 @@ FoamMappedInletBCBase::createPatchProcMap() createMapComm(getFvMesh(), face_centres, map_procs, inlet_procs, comm); if (map_procs.empty()) - mooseError("Mapped inlet invalid. No intersection of mapped plane and mesh."); + mooseError("The mapped inlet plane does not intersect the mesh."); if (_mpi_comm == MPI_COMM_NULL) // process not in mapped or inlet planes return; @@ -199,8 +227,8 @@ FoamMappedInletBCBase::createPatchProcMap() for (auto i = 0lu; i < inlet_procs.size(); ++i) { Foam::vectorField field; - Foam::UIPstream recieve(inlet_procs[i], send_points); - recieve >> field; + Foam::UIPstream receive(inlet_procs[i], send_points); + receive >> field; auto & vec = _send_map[inlet_procs[i]]; auto & recv_indices = recv_indices_procs[i]; auto indices = findIndices(field + _offset, map_comm); @@ -213,7 +241,7 @@ FoamMappedInletBCBase::createPatchProcMap() recv_indices.push_back(j); } } - if (vec.size() == 0) + if (vec.empty()) _send_map.erase(inlet_procs[i]); // Let original processes know which points will come from each rank @@ -273,7 +301,7 @@ FoamMappedInletBCBase::createPatchProcMap() } std::vector -FoamMappedInletBCBase::findIndices(const Foam::pointField & locations, const MPI_Comm & comm) +FoamMappedInletBCBase::findIndices(const Foam::pointField & locations, MPI_Comm comm) { std::vector indices(locations.size(), -1); for (auto i = 0; i < locations.size(); ++i) @@ -281,21 +309,21 @@ FoamMappedInletBCBase::findIndices(const Foam::pointField & locations, const MPI indices[i] = _mesh_searcher.findCell(locations[i], Foam::pointInCellShapes::facePlanes); } - auto bad_indices = findLocationErrors(indices, comm, FindError::NOT_FOUND); - for (auto i = 0lu; i < bad_indices.size(); ++i) + auto missing_indices = findLocatorIssues(indices, comm, LocatorIssue::NOT_FOUND); + for (auto i = 0lu; i < missing_indices.size(); ++i) { - int celli = _mesh_searcher.findNearestCell(locations[bad_indices[i]]); + int celli = _mesh_searcher.findNearestCell(locations[missing_indices[i]]); bool in_cell = Foam::pointInCell( - locations[bad_indices[i]], getFvMesh(), celli, Foam::pointInCellShapes::facePlanes); - indices[bad_indices[i]] = (in_cell) ? celli : -1; + locations[missing_indices[i]], getFvMesh(), celli, Foam::pointInCellShapes::facePlanes); + indices[missing_indices[i]] = (in_cell) ? celli : -1; } // Check not founds - std::vector notfound_points = findLocationErrors(indices, comm, FindError::NOT_FOUND); - if (notfound_points.size() > 0) + std::vector notfound_indices = findLocatorIssues(indices, comm, LocatorIssue::NOT_FOUND); + if (!notfound_indices.empty()) { std::stringstream err_msg; - for (auto idx : notfound_points) + for (auto idx : notfound_indices) { const auto & loc = locations[idx]; err_msg << "\t(" << loc.x() << ", " << loc.y() << ", " << loc.z() << ")\n"; @@ -303,64 +331,28 @@ FoamMappedInletBCBase::findIndices(const Foam::pointField & locations, const MPI mooseError("Locations not found:\n", err_msg.str()); } - // Deduplication - // check if deduplication is required - auto deconflict_points = findLocationErrors(indices, comm, FindError::MANY_FOUND); - if (deconflict_points.empty()) + // Resolve locations claimed by multiple ranks + auto conflict_indices = findLocatorIssues(indices, comm, LocatorIssue::MANY_FOUND); + if (conflict_indices.empty()) return indices; - // deconflict step 1 - // Use cell centre closest to location - std::vector dist(deconflict_points.size(), std::numeric_limits::max()); - std::vector gl_dist(deconflict_points.size(), std::numeric_limits::max()); - - for (auto i = 0lu; i < deconflict_points.size(); ++i) - { - int idx = deconflict_points[i]; - - if (indices[idx] != -1) - dist[i] = Foam::mag(locations[idx] - getFvMesh().cellCentres()[indices[idx]]); - } - - MPI_Allreduce(dist.data(), gl_dist.data(), deconflict_points.size(), MPI_DOUBLE, MPI_MIN, comm); - - for (auto i = 0lu; i < deconflict_points.size(); ++i) - { - int idx = deconflict_points[i]; - if (dist[i] != gl_dist[i]) - indices[idx] = -1; - } + // Prefer the cell whose centre is closest to the mapped location. + retainClosestCells(indices, conflict_indices, locations, getFvMesh().cellCentres(), comm); // Check again for conflicts - deconflict_points = findLocationErrors(indices, comm, FindError::MANY_FOUND); - if (deconflict_points.empty()) + conflict_indices = findLocatorIssues(indices, comm, LocatorIssue::MANY_FOUND); + if (conflict_indices.empty()) return indices; - // 2. Choose which ever is closest to the inlet - dist.assign(deconflict_points.size(), std::numeric_limits::max()); - gl_dist.assign(deconflict_points.size(), std::numeric_limits::max()); - for (auto i = 0lu; i < deconflict_points.size(); ++i) - { - int idx = deconflict_points[i]; - - if (indices[idx] != -1) - dist[i] = Foam::mag(getFvMesh().cellCentres()[indices[idx]] - (locations[idx] - _offset)); - } - - MPI_Allreduce(dist.data(), gl_dist.data(), deconflict_points.size(), MPI_DOUBLE, MPI_MIN, comm); + // Break remaining ties using the distance to the original inlet location. + const Foam::pointField inlet_locations = locations - _offset; + retainClosestCells(indices, conflict_indices, inlet_locations, getFvMesh().cellCentres(), comm); - for (auto i = 0lu; i < deconflict_points.size(); ++i) - { - int idx = deconflict_points[i]; - if (dist[i] != gl_dist[i]) - indices[idx] = -1; - } - - deconflict_points = findLocationErrors(indices, comm, FindError::MANY_FOUND); - if (!deconflict_points.empty()) + conflict_indices = findLocatorIssues(indices, comm, LocatorIssue::MANY_FOUND); + if (!conflict_indices.empty()) { std::stringstream err_msg; - for (auto idx : deconflict_points) + for (auto idx : conflict_indices) { const auto & loc = locations[idx]; err_msg << "\t(" << loc.x() << ", " << loc.y() << ", " << loc.z() << ")\n"; @@ -411,7 +403,7 @@ FoamMappedInletBCBase::getMappedArray(const Foam::word & name) Foam::PstreamBuffers sendBuf( Foam::UPstream::commsTypes::nonBlocking, Foam::UPstream::msgType(), _foam_comm); - if (_send_map.size() > 0) + if (!_send_map.empty()) { auto & var = getFvMesh().lookupObject>(name); for (auto & pair : _send_map) @@ -430,7 +422,7 @@ FoamMappedInletBCBase::getMappedArray(const Foam::word & name) sendBuf.finishedSends(true); Foam::Field boundaryData(boundary_patch.size()); - if (_recv_map.size() > 0) + if (!_recv_map.empty()) { for (auto & pair : _recv_map) { @@ -449,3 +441,30 @@ FoamMappedInletBCBase::getMappedArray(const Foam::word & name) return boundaryData; } + +Foam::label +FoamMappedInletBCBase::createCommunicator(const Foam::label parent_comm, + const std::vector & procs, + MPI_Comm & new_comm) +{ + Foam::label foam_comm; + if (Foam::UPstream::parRun()) + { + Foam::labelList foam_procs(procs.begin(), procs.end()); + foam_comm = Foam::UPstream::allocateCommunicator(parent_comm, foam_procs, true); + new_comm = Foam::PstreamGlobals::MPICommunicators_[foam_comm]; + } + else + { + foam_comm = Foam::UPstream::worldComm; + new_comm = Foam::PstreamGlobals::MPICommunicators_[Foam::UPstream::worldComm]; + } + return foam_comm; +} + +void +FoamMappedInletBCBase::destroyCommunicator(Foam::label comm) +{ + if (Foam::UPstream::parRun()) + Foam::UPstream::freeCommunicator(comm); +} diff --git a/src/bcs/FoamScalarBulkMappedInletBC.C b/src/bcs/FoamScalarBulkMappedInletBC.C index 54ca4388..555a4d27 100644 --- a/src/bcs/FoamScalarBulkMappedInletBC.C +++ b/src/bcs/FoamScalarBulkMappedInletBC.C @@ -1,14 +1,16 @@ #include "FoamScalarBulkMappedInletBC.h" #include "InputParameters.h" -#include "MooseError.h" #include "MooseTypes.h" -#include "PstreamReduceOps.H" #include "Registry.h" +#include "PstreamReduceOps.H" #include "ops.H" #include "volFieldsFwd.H" +#include #include +#include + registerMooseObject("hippoApp", FoamScalarBulkMappedInletBC); InputParameters diff --git a/test/tests/bcs/mapped_inlet/tests b/test/tests/bcs/mapped_inlet/tests index 4ef88cfa..8e18a2fe 100644 --- a/test/tests/bcs/mapped_inlet/tests +++ b/test/tests/bcs/mapped_inlet/tests @@ -12,7 +12,7 @@ min_parallel = 5 max_parallel = 5 cli_args = "FoamBCs/mass_flowx/translation_vector='-0.5 0 0'" - expect_err = "Mapped inlet invalid. No intersection of mapped plane and mesh." + expect_err = "The mapped inlet plane does not intersect the mesh." [] [run_err2] type = RunException diff --git a/test/tests/bcs/receiver_pp/test.py b/test/tests/bcs/receiver_pp/test.py index 7761bfed..0f3a6dc4 100644 --- a/test/tests/bcs/receiver_pp/test.py +++ b/test/tests/bcs/receiver_pp/test.py @@ -12,7 +12,7 @@ class TestFoamBCFixedGradient(unittest.TestCase): def test_diffusion_flux_x(self): """ Test case for imposing diffusion flux BCs using the default value of the - underlying reciever. + underlying receiver. Solves laplace equation with right BC being \partial_x T = 1. The analytical solution is x. diff --git a/test/tests/bcs/receiver_pp/tests b/test/tests/bcs/receiver_pp/tests index 5ad5468b..4f9f1179 100644 --- a/test/tests/bcs/receiver_pp/tests +++ b/test/tests/bcs/receiver_pp/tests @@ -1,6 +1,6 @@ [Tests] [receiver_pp] - requirement = "Tests the under-the-hood reciever in the Postprocessor BCs when using the default parameter" + requirement = "Tests the under-the-hood receiver in the Postprocessor BCs when using the default parameter" [setup] type = RunCommand command = 'bash -c "foamCleanCase -case foam && blockMesh -case foam"' From c8c40a74b95345771db070c9a5ebf2ed1d22808f Mon Sep 17 00:00:00 2001 From: Matthew Falcone Date: Thu, 27 Aug 2026 17:06:46 +0100 Subject: [PATCH 7/7] Abstract common code into anonymous namespace --- include/bcs/FoamMappedInletBCBase.h | 6 +- src/bcs/FoamMappedInletBCBase.C | 102 ++++++++++++++++++---------- 2 files changed, 71 insertions(+), 37 deletions(-) diff --git a/include/bcs/FoamMappedInletBCBase.h b/include/bcs/FoamMappedInletBCBase.h index 09542b1e..1b0e5523 100644 --- a/include/bcs/FoamMappedInletBCBase.h +++ b/include/bcs/FoamMappedInletBCBase.h @@ -12,7 +12,7 @@ class FoamMappedInletBCBase : public FoamPostprocessorBCBase FoamMappedInletBCBase(const InputParameters & params); - virtual ~FoamMappedInletBCBase() { destroyCommunicator(_foam_comm); } + ~FoamMappedInletBCBase() override { destroyCommunicator(_foam_comm); } protected: Foam::vector _offset; @@ -41,8 +41,8 @@ class FoamMappedInletBCBase : public FoamPostprocessorBCBase // create/assign communicators for the transfers between map and inlet planes void createMapComm(const Foam::fvMesh & mesh, const Foam::vectorField & face_centres, - std::vector & send_process, - std::vector & recv_process, + std::vector & map_processes, + std::vector & inlet_processes, MPI_Comm comm); // find indices of cells containing mapped inlet points or raise error if not found diff --git a/src/bcs/FoamMappedInletBCBase.C b/src/bcs/FoamMappedInletBCBase.C index 2c271897..ae00c900 100644 --- a/src/bcs/FoamMappedInletBCBase.C +++ b/src/bcs/FoamMappedInletBCBase.C @@ -10,11 +10,11 @@ #include #include -#include -#include #include #include #include +#include +#include namespace { @@ -79,7 +79,8 @@ findLocatorIssues(const std::vector & indices, MPI_Comm comm, LocatorIssue void retainClosestCells(std::vector & indices, const std::vector & conflict_indices, - const Foam::pointField & target_locations, + const Foam::pointField & locations, + const Foam::vector & offset, const Foam::vectorField & cell_centres, MPI_Comm comm) { @@ -90,7 +91,7 @@ retainClosestCells(std::vector & indices, { const int idx = conflict_indices[i]; if (indices[idx] >= 0) - local_distances[i] = Foam::mag(target_locations[idx] - cell_centres[indices[idx]]); + local_distances[i] = Foam::mag(locations[idx] + offset - cell_centres[indices[idx]]); } MPI_Allreduce(local_distances.data(), @@ -104,6 +105,39 @@ retainClosestCells(std::vector & indices, if (local_distances[i] != minimum_distances[i]) indices[conflict_indices[i]] = -1; } + +std::string +getLocationList(const Foam::pointField & locations, const std::vector & indices) +{ + std::stringstream msg; + for (auto idx : indices) + { + const auto & loc = locations[idx]; + msg << "\t(" << loc.x() << ", " << loc.y() << ", " << loc.z() << ")\n"; + } + return msg.str(); +} + +void +retainLowestRank(std::vector & indices, + const std::vector & conflict_indices, + MPI_Comm comm) +{ + int rank; + MPI_Comm_rank(comm, &rank); + + std::vector minimum_ranks(conflict_indices.size(), std::numeric_limits::max()); + for (auto i = 0lu; i < conflict_indices.size(); ++i) + if (indices[conflict_indices[i]] >= 0) + minimum_ranks[i] = rank; + + MPI_Allreduce( + MPI_IN_PLACE, minimum_ranks.data(), conflict_indices.size(), MPI_INT, MPI_MIN, comm); + + for (auto i = 0lu; i < conflict_indices.size(); ++i) + if (rank != minimum_ranks[i]) + indices[conflict_indices[i]] = -1; +} } bool @@ -143,8 +177,8 @@ FoamMappedInletBCBase::intersectMapPlane(const Foam::fvMesh & mesh, Real cart_bb void FoamMappedInletBCBase::createMapComm(const Foam::fvMesh & mesh, const Foam::vectorField & face_centres, - std::vector & map_process, - std::vector & inlet_process, + std::vector & map_processes, + std::vector & inlet_processes, MPI_Comm comm) { Real cart_bbox[6]; @@ -160,8 +194,8 @@ FoamMappedInletBCBase::createMapComm(const Foam::fvMesh & mesh, MPI_Allgather(&mappedPlaneProcess, 1, MPI_INT, map_procs.data(), 1, MPI_INT, comm); MPI_Allgather(&inletPlaneProcess, 1, MPI_INT, inlet_procs.data(), 1, MPI_INT, comm); - map_process.clear(); - inlet_process.clear(); + map_processes.clear(); + inlet_processes.clear(); // create list of processes in new communicator and whether they are in the inlet or mapped plane std::vector processes; @@ -172,9 +206,9 @@ FoamMappedInletBCBase::createMapComm(const Foam::fvMesh & mesh, { processes.push_back(i); if (inlet_procs[i]) - inlet_process.push_back(j); + inlet_processes.push_back(j); if (map_procs[i]) - map_process.push_back(j); + map_processes.push_back(j); ++j; } } @@ -189,7 +223,7 @@ FoamMappedInletBCBase::createPatchProcMap() auto face_centres = boundary.Cf(); std::vector map_procs, inlet_procs; - const MPI_Comm & comm{Foam::PstreamGlobals::MPICommunicators_[Foam::UPstream::worldComm]}; + const MPI_Comm comm{Foam::PstreamGlobals::MPICommunicators_[Foam::UPstream::worldComm]}; createMapComm(getFvMesh(), face_centres, map_procs, inlet_procs, comm); if (map_procs.empty()) @@ -262,14 +296,14 @@ FoamMappedInletBCBase::createPatchProcMap() if (isInletProc) // create map to determine where data from map processes should go { std::vector assignment_count(face_centres.size(), 0); - for (auto & proc : map_procs) + for (const auto & proc : map_procs) { int size{0}; MPI_Recv(&size, 1, MPI_INT, proc, 0, _mpi_comm, MPI_STATUS_IGNORE); std::vector recv_indices(size); MPI_Recv(recv_indices.data(), size, MPI_INT, proc, 1, _mpi_comm, MPI_STATUS_IGNORE); - for (auto & index : recv_indices) + for (const auto & index : recv_indices) { if (index < 0 || index >= face_centres.size()) mooseError( @@ -313,22 +347,20 @@ FoamMappedInletBCBase::findIndices(const Foam::pointField & locations, MPI_Comm for (auto i = 0lu; i < missing_indices.size(); ++i) { int celli = _mesh_searcher.findNearestCell(locations[missing_indices[i]]); + if (celli < 0) + continue; + bool in_cell = Foam::pointInCell( locations[missing_indices[i]], getFvMesh(), celli, Foam::pointInCellShapes::facePlanes); indices[missing_indices[i]] = (in_cell) ? celli : -1; } - // Check not founds + // Check whether any locations are not found std::vector notfound_indices = findLocatorIssues(indices, comm, LocatorIssue::NOT_FOUND); if (!notfound_indices.empty()) { - std::stringstream err_msg; - for (auto idx : notfound_indices) - { - const auto & loc = locations[idx]; - err_msg << "\t(" << loc.x() << ", " << loc.y() << ", " << loc.z() << ")\n"; - } - mooseError("Locations not found:\n", err_msg.str()); + auto err_msg = getLocationList(locations, notfound_indices); + mooseError("Locations not found:\n", err_msg); } // Resolve locations claimed by multiple ranks @@ -337,27 +369,29 @@ FoamMappedInletBCBase::findIndices(const Foam::pointField & locations, MPI_Comm return indices; // Prefer the cell whose centre is closest to the mapped location. - retainClosestCells(indices, conflict_indices, locations, getFvMesh().cellCentres(), comm); + retainClosestCells( + indices, conflict_indices, locations, Foam::vector::zero, getFvMesh().cellCentres(), comm); - // Check again for conflicts + // Check again for conflicts: distance criterion could leave ties conflict_indices = findLocatorIssues(indices, comm, LocatorIssue::MANY_FOUND); if (conflict_indices.empty()) return indices; // Break remaining ties using the distance to the original inlet location. - const Foam::pointField inlet_locations = locations - _offset; - retainClosestCells(indices, conflict_indices, inlet_locations, getFvMesh().cellCentres(), comm); + retainClosestCells( + indices, conflict_indices, locations, -_offset, getFvMesh().cellCentres(), comm); conflict_indices = findLocatorIssues(indices, comm, LocatorIssue::MANY_FOUND); if (!conflict_indices.empty()) { - std::stringstream err_msg; - for (auto idx : conflict_indices) - { - const auto & loc = locations[idx]; - err_msg << "\t(" << loc.x() << ", " << loc.y() << ", " << loc.z() << ")\n"; - } - mooseError("Locations found in multiple ranks\n", err_msg.str()); + int rank; + MPI_Comm_rank(comm, &rank); + if (rank == 0) + mooseWarning("Mapped locations assigned to multiple ranks after deconflicting. Using the " + "lowest rank:\n", + getLocationList(locations, conflict_indices)); + + retainLowestRank(indices, conflict_indices, comm); } return indices; @@ -406,7 +440,7 @@ FoamMappedInletBCBase::getMappedArray(const Foam::word & name) if (!_send_map.empty()) { auto & var = getFvMesh().lookupObject>(name); - for (auto & pair : _send_map) + for (const auto & pair : _send_map) { auto & proc = pair.first; auto & send_indices = pair.second; @@ -424,7 +458,7 @@ FoamMappedInletBCBase::getMappedArray(const Foam::word & name) Foam::Field boundaryData(boundary_patch.size()); if (!_recv_map.empty()) { - for (auto & pair : _recv_map) + for (const auto & pair : _recv_map) { auto & proc = pair.first; auto & recv_indices = pair.second;