Add HTC shadowing from OpenFOAM

hippo.mk

@@ -20,7 +20,9 @@ ADDITIONAL_INCLUDES := \
     -isystem $(FOAM_INCLUDE_ROOT)/OpenFOAM/lnInclude \
     -isystem $(FOAM_INCLUDE_ROOT)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -isystem $(FOAM_INCLUDE_ROOT)/ThermophysicalTransportModels/fluidThermo/lnInclude \
+    -isystem $(FOAM_INCLUDE_ROOT)/ThermophysicalTransportModels/fluid/lnInclude \
     -isystem $(FOAM_INCLUDE_ROOT)/ThermophysicalTransportModels/lnInclude \
+    -isystem $(FOAM_INCLUDE_ROOT)/ThermophysicalTransportModels/thermophysicalTransportModel/lnInclude \
     -isystem $(FOAM_INCLUDE_ROOT)/MomentumTransportModels/compressible/lnInclude \
     -isystem $(FOAM_INCLUDE_ROOT)/physicalProperties/lnInclude \
     -isystem $(FOAM_INCLUDE_ROOT)/thermophysicalModels/basic/lnInclude \

include/problems/FoamProblem.h

@@ -3,7 +3,7 @@
 #include "FoamMesh.h"
 #include "FoamPostprocessorBase.h"
 #include "FoamSolver.h"
-#include "FoamVariableField.h"
+#include "FoamFieldBase.h"
 #include "FoamBCBase.h"
 
 #include <ExternalProblem.h>
@@ -49,7 +49,7 @@ class FoamProblem : public ExternalProblem
   FoamMesh * _foam_mesh = nullptr;
   Hippo::FoamSolver _solver;
 
-  std::vector<FoamVariableField *> _foam_variables;
+  std::vector<FoamFieldBase *> _foam_variables;
   std::vector<FoamBCBase *> _foam_bcs;
   std::vector<FoamPostprocessorBase *> _foam_postprocessor;
 };

include/userobjects/AdjacentCellBulkTemperature.h

@@ -0,0 +1,33 @@
+#pragma once
+
+#include "GeneralUserObject.h"
+#include "HippoInterface.h"
+#include "InputParameters.h"
+#include "KDTree.h"
+#include <fvPatch.H>
+#include <memory>
+#include <scalarField.H>
+
+class AdjacentCellBulkTemperature : public GeneralUserObject, public HippoInterface
+{
+public:
+  static InputParameters validParams();
+  AdjacentCellBulkTemperature(const InputParameters & params);
+
+  void execute() override;
+  void initialize() override {};
+  void finalize() override {};
+
+  Real spatialValue(const Point & point) const override;
+
+private:
+  const Foam::fvPatch & getFoamPatch(const std::string & boundary);
+  void buildKDTree();
+
+  std::vector<int> _mpi_sizes;
+  std::vector<int> _mpi_displs;
+  const Foam::fvPatch & _foam_patch;
+  std::vector<Point> _kd_centres;
+  Foam::scalarField _gl_t_adjacent;
+  std::unique_ptr<KDTree> _kd_tree;
+};

include/userobjects/FixedBulkTemperature.h

@@ -0,0 +1,20 @@
+#pragma once
+
+#include "GeneralUserObject.h"
+#include "InputParameters.h"
+
+class FixedBulkTemperature : public GeneralUserObject
+{
+public:
+  static InputParameters validParams();
+  FixedBulkTemperature(const InputParameters & params);
+
+  void execute() override {}
+  void initialize() override {}
+  void finalize() override {};
+
+  Real spatialValue(const Point & point) const override;
+
+private:
+  const Real _T_bulk;
+};

include/util/hippoUtils.h

@@ -1,6 +1,9 @@
 #pragma once
 
 #include "MooseError.h"
+#include "MooseTypes.h"
+#include "MooseVariableFieldBase.h"
+#include "FoamMesh.h"
 #include <InputParameters.h>
 #include <string>
 
@@ -59,5 +62,25 @@ listFromVector(std::vector<StrType> vec, const char * sep = ", ")
   auto binary_op = [&](const std::string & acc, const std::string & it) { return acc + sep + it; };
   return std::accumulate(vec.begin() + 1, vec.end(), str, binary_op);
 }
+
+template <typename FoamField>
+inline void
+copyFieldFoamToMoose(const FoamMesh & mesh,
+                     const FoamField & field,
+                     MooseVariableFieldBase & moose_var,
+                     SubdomainID subdomain)
+{
+  size_t patch_count = mesh.getPatchCount(subdomain);
+  size_t patch_offset = mesh.getPatchOffset(subdomain);
+  for (size_t j = 0; j < patch_count; ++j)
+  {
+    auto elem = patch_offset + j;
+    auto elem_ptr = mesh.getElemPtr(elem + mesh.rank_element_offset);
+    assert(elem_ptr);
+    auto dof_t = elem_ptr->dof_number(moose_var.sys().number(), moose_var.number(), 0);
+    moose_var.sys().solution().set(dof_t, field[j]);
+  }
+}
+
 }
 }

include/variables/FoamFieldBase.h

@@ -0,0 +1,17 @@
+#pragma once
+
+#include "HippoObject.h"
+
+class FoamFieldBase : public HippoObject
+{
+public:
+  static InputParameters validParams();
+
+  explicit FoamFieldBase(const InputParameters & params);
+
+  // transfer variable from OpenFOAM field to MOOSE variable
+  virtual void transferVariable() = 0;
+
+  // returns the name of the foam variable this object shadows
+  virtual std::string foamVariable() const = 0;
+};

include/variables/FoamHeatTransferCoeff.h

@@ -0,0 +1,32 @@
+#pragma once
+
+#include "FEProblemBase.h"
+#include "FoamFieldBase.h"
+#include "InputParameters.h"
+#include "FoamMesh.h"
+#include "MooseObject.h"
+
+#include "MooseTypes.h"
+#include "fvCFD_moose.h"
+#include <fvPatchFieldsFwd.H>
+#include <scalarField.H>
+
+class FoamHeatTransferCoeff : public FoamFieldBase
+{
+public:
+  static InputParameters validParams();
+
+  explicit FoamHeatTransferCoeff(const InputParameters & params);
+
+  virtual void transferVariable();
+
+  std::string foamVariable() const { return ""; };
+
+protected:
+  Foam::scalarField calculateHTC();
+  const Foam::scalarField calculate_qw(const Foam::fvPatchScalarField & Tbf);
+  // Pointer to the FoamMesh object
+  FoamMesh & _mesh;
+  const SubdomainID _subdomain;
+  const UserObjectName & _t_bulk_uo_name;
+};

include/variables/FoamVariableField.h

@@ -1,19 +1,20 @@
 #pragma once
 
-#include "HippoObject.h"
+#include "FoamFieldBase.h"
+#include "FoamMesh.h"
 
-class FoamVariableField : public HippoObject
+class FoamVariableField : public FoamFieldBase
 {
 public:
   static InputParameters validParams();
 
   explicit FoamVariableField(const InputParameters & params);
 
   // transfer variable from OpenFOAM field to MOOSE variable
-  virtual void transferVariable();
+  virtual void transferVariable() override;
 
   // returns the name of the foam variable this object shadows
-  std::string foamVariable() const { return _foam_variable; };
+  virtual std::string foamVariable() const override { return _foam_variable; };
 
 protected:
   // variable name or functionObject to be shadowed

src/actions/AddFoamVariableAction.C

@@ -1,6 +1,5 @@
 #include "AddFoamVariableAction.h"
-#include "FoamProblem.h"
-#include "FoamVariableField.h"
+#include "FoamFieldBase.h"
 #include "hippoUtils.h"
 
 #include "InputParameters.h"
@@ -31,7 +30,7 @@ AddFoamVariableAction::act()
   {
     createAuxVariable();
 
-    _problem->addObject<FoamVariableField>(_type, _name, _moose_object_pars, false);
+    _problem->addObject<FoamFieldBase>(_type, _name, _moose_object_pars, false);
   }
 }
 

src/problems/FoamProblem.C

@@ -1,14 +1,14 @@
-#include "FoamVariableField.h"
 #include "FoamMesh.h"
 #include "FoamProblem.h"
 #include "FoamSolver.h"
 #include "hippoUtils.h"
 #include "FoamSidePostprocessor.h"
 
-#include "Attributes.h"
-#include "ExternalProblem.h"
-#include "VariadicTable.h"
-#include "MooseTypes.h"
+#include <AuxiliarySystem.h>
+#include <MooseError.h>
+#include <MooseTypes.h>
+#include <MooseVariableFieldBase.h>
+#include <algorithm>
 #include "InputParameters.h"
 #include "VariadicTable.h"
 

src/userobjects/AdjacentCellBulkTemperature.C

@@ -0,0 +1,144 @@
+#include "AdjacentCellBulkTemperature.h"
+#include "FoamProblem.h"
+#include "MooseTypes.h"
+#include "petsclog.h"
+#include <Field.H>
+#include <ListOps.H>
+#include <Pstream/mpi/PstreamGlobals.H>
+#include <PstreamReduceOps.H>
+#include <UList.H>
+#include <UPstream.H>
+#include <VectorSpace.H>
+#include <face.H>
+#include <fvMesh.H>
+#include <memory>
+#include <ops.H>
+#include <scalar.H>
+#include <scalarField.H>
+#include <vector>
+#include <vectorField.H>
+#include <KDTree.h>
+
+registerMooseObject("hippoApp", AdjacentCellBulkTemperature);
+
+InputParameters
+AdjacentCellBulkTemperature::validParams()
+{
+  InputParameters params = GeneralUserObject::validParams();
+  params.addRequiredParam<SubdomainName>("boundary", "Boundary this object applies to.");
+  params.addParam<std::string>("T_name", "T", "Temperature name in the OpenFOAM simulation.");
+  params.addParam<int>("max_leaf_size", 10, "Maximum leaf size to use for the KD tree search.");
+  params.addParam<bool>(
+      "reconstruct_tree", false, "Whether to rebuild the KD tree on each execution.");
+  return params;
+}
+
+AdjacentCellBulkTemperature::AdjacentCellBulkTemperature(const InputParameters & params)
+  : GeneralUserObject(params),
+    HippoInterface(this),
+    _foam_patch{getFoamPatch(getParam<SubdomainName>("boundary"))},
+    _kd_centres{},
+    _gl_t_adjacent{},
+    _kd_tree{}
+
+{
+}
+
+void
+AdjacentCellBulkTemperature::buildKDTree()
+{
+  Foam::vectorField face_centres = _foam_patch.Cf();
+  const int l_num_faces = face_centres.size();
+  const int nProcs{Foam::UPstream::nProcs()};
+  const MPI_Comm comm{(nProcs == 1) ? MPI_COMM_WORLD : Foam::PstreamGlobals::MPI_COMM_FOAM};
+
+  if (_gl_t_adjacent.size() == 0)
+  {
+    int gl_size;
+    MPI_Allreduce(&l_num_faces, &gl_size, 1, MPI_INT, MPI_SUM, comm);
+    _gl_t_adjacent.resize(gl_size);
+    _kd_centres.resize(gl_size);
+  }
+
+  std::vector<int> face_sizes(nProcs);
+  std::vector<int> face_displs(nProcs);
+  _mpi_sizes.resize(nProcs);
+  _mpi_displs.resize(nProcs);
+
+  // Create parameters for MPI gather operations to be reused in execute
+  MPI_Allgather(&l_num_faces, 1, MPI_INT, _mpi_sizes.data(), 1, MPI_INT, comm);
+  _mpi_displs[0] = 0;
+  for (int i = 0; i < nProcs - 1; ++i)
+  {
+    _mpi_displs[i + 1] = _mpi_displs[i] + _mpi_sizes[i];
+  }
+
+  // Account for vector for gathering coordinates
+  for (int i = 0; i < nProcs; ++i)
+  {
+    face_sizes[i] = _mpi_sizes[i] * 3;
+    face_displs[i] = _mpi_displs[i] * 3;
+  }
+
+  Foam::vectorField gl_face_centres(_gl_t_adjacent.size());
+  MPI_Allgatherv(face_centres.data(),
+                 face_centres.size() * 3,
+                 MPI_DOUBLE,
+                 gl_face_centres.data(),
+                 face_sizes.data(),
+                 face_displs.data(),
+                 MPI_DOUBLE,
+                 comm);
+
+  // Create KD tree
+  for (int i = 0; i < gl_face_centres.size(); ++i)
+  {
+    const Foam::point & p = gl_face_centres[i];
+    _kd_centres[i] = Point(p.x(), p.y(), p.z());
+  }
+  _kd_tree = std::make_unique<KDTree>(_kd_centres, getParam<int>("max_leaf_size"));
+}
+
+void
+AdjacentCellBulkTemperature::execute()
+{
+  if (getParam<bool>("reconstruct_tree") || _gl_t_adjacent.size() == 0)
+    buildKDTree();
+
+  // Gather Temperature on each execute
+  Foam::scalarField l_t_adjacent{_foam_patch.Cf().size()};
+
+  const auto & T = getFvMesh().lookupObject<Foam::volScalarField>(getParam<std::string>("T_name"));
+  for (int i = 0; i < l_t_adjacent.size(); ++i)
+  {
+    int idx = _foam_patch.faceCells()[i];
+    l_t_adjacent[i] = T[idx];
+  }
+
+  const int nProcs{Foam::UPstream::nProcs()};
+  const MPI_Comm comm{(nProcs == 1) ? MPI_COMM_WORLD : Foam::PstreamGlobals::MPI_COMM_FOAM};
+  MPI_Allgatherv(l_t_adjacent.data(),
+                 l_t_adjacent.size(),
+                 MPI_DOUBLE,
+                 _gl_t_adjacent.data(),
+                 _mpi_sizes.data(),
+                 _mpi_displs.data(),
+                 MPI_DOUBLE,
+                 comm);
+}
+
+const Foam::fvPatch &
+AdjacentCellBulkTemperature::getFoamPatch(const std::string & boundary)
+{
+  const Foam::label & boundaryId{getFvMesh().boundary().findIndex(boundary)};
+  return getFvMesh().boundary()[boundaryId];
+}
+
+Real
+AdjacentCellBulkTemperature::spatialValue(const Point & point) const
+{
+  const Foam::point p_moose{point(0), point(1), point(2)};
+  std::vector<size_t> idx_vec;
+  _kd_tree->neighborSearch(point, 1, idx_vec);
+  return _gl_t_adjacent[idx_vec.at(0)];
+}