Interface for Python 1d flow code to axisymmetric EM solver

.gitattributes

@@ -27,3 +27,4 @@ test/ref_solns/interpInlet/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -
 test/meshes/spongeBox.msh filter=lfs diff=lfs merge=lfs -text
 test/ref_solns/sgsLoMach/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -text
 test/ref_solns/aveLoMach/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -text
+test/ref_solns/flow1d_coupling.h5 filter=lfs diff=lfs merge=lfs -text

docker/test/Dockerfile

@@ -30,6 +30,7 @@ RUN yum -y install boost-gnu9-mpich-ohpc
 
 RUN yum -y install python38-pip
 RUN pip3 install matplotlib
+RUN pip3 install h5py
 
 RUN yum -y install procps-ng
 RUN yum -y install diffutils
@@ -150,7 +151,7 @@ RUN yum -y install gsl-gnu9-ohpc
 RUN pip3 install "pybind11[global]"
 
 RUN . /etc/profile.d/lmod.sh \
-    && pip3 install mpi4py
+    && pip3 install mpi4py==3.1.5
 
 # Register new libs installed into /usr/local/lib with linker
 RUN echo "/usr/local/lib" > /etc/ld.so.conf.d/class.conf

src/Makefile.am

@@ -36,8 +36,8 @@ cxx_sources          = averaging.cpp faceGradientIntegration.cpp M2ulPhyS.cpp rh
 		       dgNonlinearForm.cpp gradients.cpp gradNonLinearForm.cpp quasimagnetostatic.cpp tps.cpp chemistry.cpp \
                        reaction.cpp collision_integrals.cpp argon_transport.cpp source_term.cpp gpu_constructor.cpp \
                        independent_coupling.cpp cycle_avg_joule_coupling.cpp table.cpp radiation.cpp lte_mixture.cpp \
-                       lte_transport_properties.cpp mixing_length_transport.cpp tps2Boltzmann.cpp  M2ulPhyS2Boltzmann.cpp \
-                       pybindings.cpp gslib_interpolator.cpp loMach.cpp loMachIO.cpp loMach_options.cpp calorically_perfect.cpp \
+                       lte_transport_properties.cpp mixing_length_transport.cpp tps2Boltzmann.cpp  M2ulPhyS2Boltzmann.cpp qms2flow1d.cpp \
+                    data_exchange_utils.cpp pybindings.cpp gslib_interpolator.cpp loMach.cpp loMachIO.cpp loMach_options.cpp calorically_perfect.cpp \
 	               thermo_chem_base.cpp split_flow_base.cpp tomboulides.cpp dirichlet_bc_helper.cpp \
                        sponge_base.cpp geometricSponge.cpp gaussianInterpExtData.cpp turb_model_base.cpp \
                        algebraicSubgridModels.cpp algebraic_rans.cpp mesh_base.cpp lte_thermo_chem.cpp cases.cpp \
@@ -53,8 +53,8 @@ headers              = averaging.hpp equation_of_state.hpp transport_properties.
 		       quasimagnetostatic.hpp gradients.hpp logger.hpp solver.hpp tps.hpp chemistry.hpp reaction.hpp \
 	               collision_integrals.hpp argon_transport.hpp source_term.hpp tps_mfem_wrap.hpp gpu_constructor.hpp \
                        independent_coupling.hpp cycle_avg_joule_coupling.hpp table.hpp radiation.hpp lte_mixture.hpp \
-                       lte_transport_properties.hpp mixing_length_transport.hpp tps2Boltzmann.hpp gslib_interpolator.hpp \
-                       loMach.hpp loMach_options.hpp calorically_perfect.hpp \
+                       lte_transport_properties.hpp mixing_length_transport.hpp data_exchange_utils.hpp tps2Boltzmann.hpp gslib_interpolator.hpp \
+                       qms2flow1d.hpp loMach.hpp loMach_options.hpp calorically_perfect.hpp \
 		       thermo_chem_base.hpp split_flow_base.hpp tomboulides.hpp dirichlet_bc_helper.hpp \
                        sponge_base.hpp geometricSponge.hpp externalData_base.hpp gaussianInterpExtData.hpp \
                        turb_model_base.hpp algebraicSubgridModels.hpp algebraic_rans.hpp mesh_base.hpp \

src/data_exchange_utils.cpp

@@ -0,0 +1,68 @@
+// -----------------------------------------------------------------------------------bl-
+// BSD 3-Clause License
+//
+// Copyright (c) 2020-2022, The PECOS Development Team, University of Texas at Austin
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// 1. Redistributions of source code must retain the above copyright notice, this
+//    list of conditions and the following disclaimer.
+//
+// 2. Redistributions in binary form must reproduce the above copyright notice,
+//    this list of conditions and the following disclaimer in the documentation
+//    and/or other materials provided with the distribution.
+//
+// 3. Neither the name of the copyright holder nor the names of its
+//    contributors may be used to endorse or promote products derived from
+//    this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+// -----------------------------------------------------------------------------------el-
+
+#include "data_exchange_utils.hpp"
+
+#ifdef HAVE_PYTHON
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+namespace tps_wrappers {
+void data_exchange_utils(py::module &m) {
+  // Can by read in numpy as np.array(data_instance, copy = False)
+  py::class_<TPS::CPUDataRead>(m, "CPUDataRead", py::buffer_protocol())
+      .def_buffer([](TPS::CPUDataRead &d) -> py::buffer_info {
+        return py::buffer_info(d.data(),                                /*pointer to buffer*/
+                               sizeof(double),                          /*size of one element*/
+                               py::format_descriptor<double>::format(), /*python struct-style format descriptor*/
+                               1,                                       /*number of dimensions*/
+                               {d.size()},                              /*buffer dimension(s)*/
+                               {d.stride() * sizeof(const double)},     /*Stride in bytes for each index*/
+                               true /*read only*/);
+      });
+
+  py::class_<TPS::CPUData>(m, "CPUData", py::buffer_protocol()).def_buffer([](TPS::CPUData &d) -> py::buffer_info {
+    return py::buffer_info(d.data(),                                /*pointer to buffer*/
+                           sizeof(double),                          /*size of one element*/
+                           py::format_descriptor<double>::format(), /*python struct-style format descriptor*/
+                           1,                                       /*number of dimensions*/
+                           {d.size()},                              /*buffer dimension(s)*/
+                           {d.stride() * sizeof(const double)},     /*Stride in bytes for each index*/
+                           false /*writetable*/);
+  });
+}
+}  // namespace tps_wrappers
+
+#endif

src/data_exchange_utils.hpp

@@ -0,0 +1,67 @@
+// -----------------------------------------------------------------------------------bl-
+// BSD 3-Clause License
+//
+// Copyright (c) 2020-2022, The PECOS Development Team, University of Texas at Austin
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// 1. Redistributions of source code must retain the above copyright notice, this
+//    list of conditions and the following disclaimer.
+//
+// 2. Redistributions in binary form must reproduce the above copyright notice,
+//    this list of conditions and the following disclaimer in the documentation
+//    and/or other materials provided with the distribution.
+//
+// 3. Neither the name of the copyright holder nor the names of its
+//    contributors may be used to endorse or promote products derived from
+//    this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+// -----------------------------------------------------------------------------------el-
+
+#include "tps.hpp"
+
+#ifndef DATA_EXCHANGE_UTILS_HPP_
+#define DATA_EXCHANGE_UTILS_HPP_
+
+namespace TPS {
+
+class CPUDataRead {
+ public:
+  CPUDataRead(const mfem::Vector &v) : data_(v.HostRead()), size_(v.Size()), stride_(1) {}
+  double *data() const { return const_cast<double *>(data_); }
+  size_t size() const { return size_; }
+  size_t stride() const { return stride_; }
+
+ private:
+  const double *data_;
+  size_t size_;
+  size_t stride_;
+};
+
+class CPUData {
+ public:
+  CPUData(mfem::Vector &v, bool rw) : data_(rw ? v.HostReadWrite() : v.HostWrite()), size_(v.Size()), stride_(1) {}
+  double *data() { return data_; }
+  size_t size() const { return size_; }
+  size_t stride() const { return stride_; }
+
+ private:
+  double *data_;
+  size_t size_;
+  size_t stride_;
+};
+}  // namespace TPS
+
+#endif  // DATA_EXCHANGE_UTILS_HPP_

src/gslib_interpolator.cpp

@@ -46,6 +46,8 @@ InterpolatorBase::InterpolatorBase() {
 
 InterpolatorBase::~InterpolatorBase() {
 #ifdef HAVE_GSLIB
+  //  Frees internal data used by FindPointsGSLIB
+  finder_->FreeData();
   delete finder_;
 #endif
 }
@@ -204,3 +206,54 @@ void PlaneInterpolator::writeAscii(std::string oname, bool rank0) const {
   // Rest of the write is handled by the base class
   InterpolatorBase::writeAscii(oname, rank0);
 }
+
+LineInterpolator::LineInterpolator(int n) : InterpolatorBase(), n_(n) {
+  start_.SetSize(0);
+  end_.SetSize(0);
+}
+
+LineInterpolator::LineInterpolator(mfem::Vector start, mfem::Vector end, int n)
+    : InterpolatorBase(), start_(start), end_(end), n_(n) {}
+
+LineInterpolator::~LineInterpolator() {}
+
+void LineInterpolator::initializeFinder(mfem::ParMesh *mesh) {
+  InterpolatorBase::initializeFinder(mesh);
+  start_.SetSize(dim_);
+  end_.SetSize(dim_);
+}
+
+void LineInterpolator::setInterpolationPoints() {
+  xyz_.SetSize(n_*dim_);
+
+  //  Vector step along line
+  mfem::Vector step = end_;
+  step -= start_;
+  step /= n_ - 1;
+
+  mfem::Vector point = start_;
+
+  for (int i = 0; i < n_; i++) {
+    xyz_[i + 0*n_] = point[0];
+    xyz_[i + 1*n_] = point[1];
+    if (dim_ == 3) {
+      xyz_[i + 2*n_] = point[2];
+    }
+    point += step;
+  }
+}
+
+void LineInterpolator::writeAscii(std::string oname, bool rank0) const {
+  if (rank0) {
+    std::cout << " Writing line data to " << oname << endl;
+
+    // First write a little info about the line
+    std::ofstream outfile;
+    outfile.open(oname, std::ios_base::app);
+    outfile << "#line start " << start_[0] << " " << start_[1] << " " << start_[2] << endl;
+    outfile << "#line end " << end_[0] << " " << end_[1] << " " << end_[2] << endl;
+    outfile.close();
+  }
+  // Rest of the write is handled by the base class
+  InterpolatorBase::writeAscii(oname, rank0);
+}

src/gslib_interpolator.hpp

@@ -86,6 +86,9 @@ class InterpolatorBase {
 
   /** Write the data to an ascii file */
   virtual void writeAscii(std::string oname, bool rank0 = true) const;
+
+  mfem::Vector &getSolution() { return soln_; }
+  mfem::Vector &getInterpolationPoints() { return xyz_; }
 };
 
 class PlaneInterpolator : public InterpolatorBase {
@@ -110,4 +113,35 @@ class PlaneInterpolator : public InterpolatorBase {
   void writeAscii(std::string oname, bool rank0 = true) const final;
 };
 
+class LineInterpolator : public InterpolatorBase {
+ private:
+  //  Endpoints of line
+  mfem::Vector start_;
+  mfem::Vector end_;
+
+  //  Number of points on line
+  int n_;
+
+ public:
+  LineInterpolator(int n = 10);
+  LineInterpolator(mfem::Vector start, mfem::Vector end, int n);
+  virtual ~LineInterpolator();
+
+  void initializeFinder(mfem::ParMesh *mesh);
+
+  void setn(int n) { n_ = n; }
+  void setStart(mfem::Vector start) { start_ = start; }
+  mfem::Vector *getStart() { return &start_; }
+  void setEnd(mfem::Vector end) { end_ = end; }
+  mfem::Vector *getEnd() { return &end_; }
+
+  /** Sets evenly spaced interpolation points
+   * 
+   * point_i = start_ + (end_ - start_)*i/(n_ - 1)
+   */
+  void setInterpolationPoints() final;
+
+  void writeAscii(std::string oname, bool rank0 = true) const final;
+};
+
 #endif  // GSLIB_INTERPOLATOR_HPP_

src/pybindings.cpp

@@ -43,13 +43,17 @@ namespace py = pybind11;
 namespace tps_wrappers {
 void tps(py::module& m);
 void tps2bolzmann(py::module& m);
+void data_exchange_utils(py::module& m);
+void qms2flow1d(py::module& m);
 }  // namespace tps_wrappers
 
 PYBIND11_MODULE(libtps, m) {
   m.doc() = "TPS Python Interface";
 
   tps_wrappers::tps(m);
   tps_wrappers::tps2bolzmann(m);
+  tps_wrappers::data_exchange_utils(m);
+  tps_wrappers::qms2flow1d(m);
 }
 
 #endif