diff --git a/.gitattributes b/.gitattributes index dda8bcaaf..e99208be6 100644 --- a/.gitattributes +++ b/.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 diff --git a/docker/test/Dockerfile b/docker/test/Dockerfile index d3869b629..50e665964 100644 --- a/docker/test/Dockerfile +++ b/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 diff --git a/src/Makefile.am b/src/Makefile.am index 423a717d5..eee7cf4cb 100644 --- a/src/Makefile.am +++ b/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 \ diff --git a/src/data_exchange_utils.cpp b/src/data_exchange_utils.cpp new file mode 100644 index 000000000..f378fb015 --- /dev/null +++ b/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 +#include +#include + +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_(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::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_(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::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 diff --git a/src/data_exchange_utils.hpp b/src/data_exchange_utils.hpp new file mode 100644 index 000000000..cb88c2195 --- /dev/null +++ b/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(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_ diff --git a/src/gslib_interpolator.cpp b/src/gslib_interpolator.cpp index 5046e57e6..1408c75f6 100644 --- a/src/gslib_interpolator.cpp +++ b/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); +} diff --git a/src/gslib_interpolator.hpp b/src/gslib_interpolator.hpp index 44bb1ae81..d7d81364b 100644 --- a/src/gslib_interpolator.hpp +++ b/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_ diff --git a/src/pybindings.cpp b/src/pybindings.cpp index 2377b127c..074d0821d 100644 --- a/src/pybindings.cpp +++ b/src/pybindings.cpp @@ -43,6 +43,8 @@ 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) { @@ -50,6 +52,8 @@ PYBIND11_MODULE(libtps, m) { tps_wrappers::tps(m); tps_wrappers::tps2bolzmann(m); + tps_wrappers::data_exchange_utils(m); + tps_wrappers::qms2flow1d(m); } #endif diff --git a/src/qms2flow1d.cpp b/src/qms2flow1d.cpp new file mode 100644 index 000000000..f1d36e07e --- /dev/null +++ b/src/qms2flow1d.cpp @@ -0,0 +1,289 @@ +// -----------------------------------------------------------------------------------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 "qms2flow1d.hpp" +#include "data_exchange_utils.hpp" + +#ifdef HAVE_PYTHON +#include +#include +#include +#endif + +#include +#include +#include + +namespace TPS { + +Qms2Flow1d::Qms2Flow1d(Tps *tps) { + tpsP_ = tps; + em_opts = nullptr; + qmsa = nullptr; + l_interp = nullptr; + + n_1d = 0; + cond_1d = nullptr; + joule_1d = nullptr; + z_coords_1d = nullptr; + radius_1d = nullptr; + n_interp = 0; +} + +Qms2Flow1d::~Qms2Flow1d() { + delete cond_1d; + delete joule_1d; + delete z_coords_1d; + delete radius_1d; + + delete em_opts; + delete qmsa; + delete l_interp; +} + +void Qms2Flow1d::initialize(const int n_1d_) { + n_1d = n_1d_; + cond_1d = new Vector(n_1d); + joule_1d = new Vector(n_1d); + z_coords_1d = new Vector(n_1d); + radius_1d = new Vector(n_1d); + + em_opts = new ElectromagneticOptions(); + qmsa = new QuasiMagnetostaticSolverAxiSym(*em_opts, tpsP_); + l_interp = new LineInterpolator(); + + n_interp = 100; + + qmsa->parseSolverOptions(); + qmsa->initialize(); + // Ensuring nodes is required for GSLib interpolation + qmsa->getMesh()->EnsureNodes(); +} + +void Qms2Flow1d::solve() { + set_plasma_conductivity_2d(); + qmsa->solve(); + set_joule_heating_1d(); +} + +void Qms2Flow1d::print_all_1d() { + std::cout << "\n--- 1d Coordinates ---" << endl; + z_coords_1d->Print(); + std::cout << "\n--- Torch Radius ---" << endl; + radius_1d->Print(); + std::cout << "\n--- 1d Plasma Conductivity ---" << endl; + cond_1d->Print(); + std::cout << "\n--- 1d Joule Heating ---" << endl; + joule_1d->Print(); +} + +double Qms2Flow1d::radial_profile(double r, double torch_r, double r_c) const { + // Constant over [0, r_c]. Zero beyond torch_r + // Cosine for continuous connection between the two + assert(0 <= r && 0 <= torch_r && 0 <= r_c); + double out = 0; + if (0 <= r && r <= r_c) { + out = 1; + } else if (r_c < r && r < torch_r) { + out = 1 + cos(M_PI*(r - r_c)/(torch_r - r_c)); + out *= 0.5; + } + return out; +} + +int Qms2Flow1d::find_z_interval(double z) const { + int i_min; + if (z < z_coords_1d->Elem(1)) { + i_min = 0; + } else if (z > z_coords_1d->Elem(n_1d - 1)) { + i_min = n_1d - 2; + } else { + i_min = 1; + int i_max = n_1d - 2; + while (i_max > i_min + 1) { + int i_test = round(0.5*(i_max - i_min) + i_min); + if (z == z_coords_1d->Elem(i_test)) { + i_min = i_test; + break; + } else if (z < z_coords_1d->Elem(i_test)) { + i_max = i_test; + } else { + i_min = i_test; + } + } + } + return i_min; +} + +double Qms2Flow1d::interpolate_z(double z, const Vector &values) const { + const int ind = find_z_interval(z); + double interp = (values[ind+1] - values[ind])/(z_coords_1d->Elem(ind+1) - z_coords_1d->Elem(ind)); + interp *= (z - z_coords_1d->Elem(ind)); + interp += values[ind]; + return interp; +} + +double Qms2Flow1d::expand_cond_2d(double r, double z) const { + const double sigma_1d = interpolate_z(z, *cond_1d); + const double torch_radius = interpolate_z(z, *radius_1d); + double sigma = radial_profile(r, torch_radius, 0.7*torch_radius); + sigma *= sigma_1d; + return sigma; +} + +void Qms2Flow1d::set_plasma_conductivity_2d() { + ParFiniteElementSpace *fes = qmsa->getFESpace(); + ParMesh *pmesh = qmsa->getMesh(); + + ParFiniteElementSpace *coord_fes = new ParFiniteElementSpace(pmesh, fes->FEColl(), 2); + ParGridFunction *coordsDof = new ParGridFunction(coord_fes); + pmesh->GetNodes(*coordsDof); + + ParGridFunction *pc = qmsa->getPlasmaConductivityGF(); + double *plasma_conductivity_gf = pc->HostWrite(); + + const int n_nodes = fes->GetNDofs(); + double r, z; + for (int i = 0; i < n_nodes; i++) { + r = coordsDof->Elem(i); + z = coordsDof->Elem(i + n_nodes); + plasma_conductivity_gf[i] = expand_cond_2d(r, z); + } + + delete coord_fes; + delete coordsDof; +} + +double Qms2Flow1d::linear_integral(Vector values, Vector coords) { + assert(values.Size() == coords.Size()); + const int n_minus_one = values.Size() - 1; + + // x_i - x_i-1 + Vector temp = Vector(coords, 1, n_minus_one); + Vector diff = temp; + temp = Vector(coords, 0, n_minus_one); + diff -= temp; + + // (v_i + v_i-1)/2 * (x_i - x_i-1) + temp = Vector(values, 0, n_minus_one); + Vector summands = temp; + temp = Vector(values, 1, n_minus_one); + summands += temp; + summands *= 0.5; + summands *= diff; + + return summands.Sum(); +} + +// Includes 2pi r Jacobian +double Qms2Flow1d::radial_integral(Vector values, Vector r_coords) { + Vector input = values; + input *= r_coords; + return 2*M_PI*linear_integral(input, r_coords); +} + +void Qms2Flow1d::set_joule_heating_1d() { + double z, torch_r; + l_interp->setn(n_interp); + l_interp->initializeFinder(qmsa->getMesh()); + + for (int i = 0; i < n_1d; i++) { + z = z_coords_1d->Elem(i); + + if (z > length) { + joule_1d->Elem(i) = 0; + } else { + torch_r = radius_1d->Elem(i); + + l_interp->getStart()->Elem(0) = 0; + l_interp->getStart()->Elem(1) = z; + l_interp->getEnd()->Elem(0) = torch_r; + l_interp->getEnd()->Elem(1) = z; + + l_interp->setInterpolationPoints(); + l_interp->interpolate(qmsa->getJouleHeatingGF()); + + Vector values = l_interp->getSolution(); + Vector r_coords = Vector(l_interp->getInterpolationPoints(), 0, n_interp); + + joule_1d->Elem(i) = radial_integral(values, r_coords); + } + } +} + +} // namespace TPS + +#ifdef HAVE_PYTHON +namespace py = pybind11; + +namespace tps_wrappers { +void qms2flow1d(py::module &m) { + py::class_(m, "Qms2Flow1d") + .def(py::init()) + .def("initialize", &TPS::Qms2Flow1d::initialize, "Initialize the 1d flow-axisymmetric EM interface", + py::arg("n_1d")) + .def("solve", &TPS::Qms2Flow1d::solve, "Solve the 1d Joule heating based on the 1d plasma conductivity") + .def("print_all", &TPS::Qms2Flow1d::print_all_1d, + "Print the 1d z coordinates, torch radius, plasma conductivity, and Joule heating from C++") + .def("total_joule_2d", &TPS::Qms2Flow1d::total_joule_heating_2d, + "Calculate the total Joule heating from the 2d solution") + .def("total_joule_1d", &TPS::Qms2Flow1d::total_joule_heating_1d, + "Calculate the total Joule heating from the 1d solution") + .def("calculate_joule", &TPS::Qms2Flow1d::set_joule_heating_1d, + "Calculates the 1d Joule heating without resolving the EM problem") + .def("set_n_interp", &TPS::Qms2Flow1d::set_n_interp, + "Sets number of interpolation points for radial integrals", + py::arg("n_interp")) + .def("set_length", &TPS::Qms2Flow1d::set_length, + "Set length of plasma torch (Joule heating is zero in the plume)", + py::arg("length")) + .def("Coordinates1d", + [](TPS::Qms2Flow1d &interface) { + return std::unique_ptr(new TPS::CPUData(interface.Coordinates1d(), false)); + }, "1d z coordinates memory location buffered read/write access") + .def("TorchRadius1d", + [](TPS::Qms2Flow1d &interface) { + return std::unique_ptr(new TPS::CPUData(interface.TorchRadius1d(), false)); + }, "1d torch radius memory location buffered read/write access") + .def("PlasmaConductivity1d", + [](TPS::Qms2Flow1d &interface) { + return std::unique_ptr(new TPS::CPUData(interface.PlasmaConductivity1d(), false)); + }, "1d plasma conductivity memory location buffered read/write access") + .def("JouleHeating1d", + [](TPS::Qms2Flow1d &interface) { + return std::unique_ptr(new TPS::CPUDataRead(interface.JouleHeating1d())); + }, "1d Joule heating memory location buffered read access"); +} + +} // namespace tps_wrappers + +#endif diff --git a/src/qms2flow1d.hpp b/src/qms2flow1d.hpp new file mode 100644 index 000000000..489878fc2 --- /dev/null +++ b/src/qms2flow1d.hpp @@ -0,0 +1,117 @@ +// -----------------------------------------------------------------------------------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- + +#ifndef QMS2FLOW1D_HPP_ +#define QMS2FLOW1D_HPP_ + +#include "tps_mfem_wrap.hpp" +#include "tps.hpp" +#include "em_options.hpp" +#include "quasimagnetostatic.hpp" +#include "gslib_interpolator.hpp" + +namespace TPS { + +class Qms2Flow1d { + public: + Qms2Flow1d(Tps *tps); + ~Qms2Flow1d(); + // Creates vectors of required size for 1d flow and initializes EM solver + void initialize(const int n_1d_); + // Expands 1d conductivity, solves EM, and reduces 2d Joule heating + void solve(); + // Prints the four 1d vectors + void print_all_1d(); + + const Vector &PlasmaConductivity1d() const { return *cond_1d; } + Vector &PlasmaConductivity1d() { return *cond_1d; } + + const Vector &JouleHeating1d() const { return *joule_1d; } + Vector &JouleHeating1d() { return *joule_1d; } + + const Vector &Coordinates1d() const { return *z_coords_1d; } + Vector &Coordinates1d() { return *z_coords_1d; } + + const Vector &TorchRadius1d() const { return *radius_1d; } + Vector &TorchRadius1d() { return *radius_1d; } + + void set_n_interp(const int n_interp_) { n_interp = n_interp_; } + void set_length(const double length_) { length = length_; } + + // Reduces 2d Joule heating by taking radial totals at 1d coordinates + void set_joule_heating_1d(); + + // Useful for comparison/convergence + double total_joule_heating_1d() { return linear_integral(*joule_1d, *z_coords_1d); } + double total_joule_heating_2d() { return qmsa->totalJouleHeating(); } + + private: + Tps *tpsP_; + ElectromagneticOptions *em_opts; + QuasiMagnetostaticSolverAxiSym *qmsa; + // Number of points in 1d flow solver + int n_1d; + Vector *cond_1d; + Vector *joule_1d; + // 1d coordinates are assumed to be sorted + Vector *z_coords_1d; + // Torch radius at z coordinates + Vector *radius_1d; + // Plasma torch length [m]. Joule heating is zero for z > L + double length = 0.315; + + // GSLib based interpolator from FE solution + LineInterpolator *l_interp; + // Number of interpolation points for radial integrals + int n_interp; + // Trapezoidal rule + double linear_integral(Vector values, Vector coords); + double radial_integral(Vector values, Vector r_coords); + + // r: radial location, torch_r: torch radius + // r_c: modeling parameter 0 <= r_c <= R + double radial_profile(double r, double torch_r, double r_c) const; + // Binary search to find interval + int find_z_interval(double z) const; + // Interpolation of 1d values along centerline via binary search + // Unrelated to GSLib interpolation of FE solution + double interpolate_z(double z, const Vector &values) const; + // Multiplies interpolant by radial profile + double expand_cond_2d(double r, double z) const; + + // Expands 1d conductivity into 2d via linear interpolation + // along z-axis multiplied by a radial profile + void set_plasma_conductivity_2d(); +}; + +} // namespace TPS +#endif // QMS2FLOW1D_HPP_ diff --git a/src/tps2Boltzmann.cpp b/src/tps2Boltzmann.cpp index 689fded06..b1e2b2771 100644 --- a/src/tps2Boltzmann.cpp +++ b/src/tps2Boltzmann.cpp @@ -34,6 +34,7 @@ */ #include "tps2Boltzmann.hpp" +#include "data_exchange_utils.hpp" #ifdef HAVE_PYTHON #include @@ -48,32 +49,6 @@ namespace TPS { -class CPUDataRead { - public: - CPUDataRead(const mfem::Vector &v) : data_(v.HostRead()), size_(v.Size()), stride_(1) {} - double *data() const { return const_cast(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_; -}; - void idenity_fun(const Vector &x, Vector &out) { for (int i(0); i < x.Size(); ++i) out[i] = x[i]; } @@ -292,28 +267,6 @@ namespace py = pybind11; namespace tps_wrappers { void tps2bolzmann(py::module &m) { - // Can by read in numpy as np.array(data_instance, copy = False) - py::class_(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::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_(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::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*/); - }); - py::enum_(m, "t2bIndex") .value("SpeciesDensities", TPS::Tps2Boltzmann::Index::SpeciesDensities) .value("ElectricField", TPS::Tps2Boltzmann::Index::ElectricField) diff --git a/test/Makefile.am b/test/Makefile.am index 0366d3e99..9ee1c876e 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -35,6 +35,7 @@ EXTRA_DIST = tap-driver.sh test_tps_splitcomm.py soln_differ inputs meshes lte- ref_solns/reactBinDiff/*.h5 \ ref_solns/reactSingleRx/*.h5 \ ref_solns/reactTable/*.h5 \ + ref_solns/flow1d_coupling.h5 \ vpath.sh die.sh count_gpus.sh sniff_mpirun.sh \ cyl3d.gpu.test cyl3d.mflow.gpu.test wedge.gpu.test \ averaging.gpu.test cyl3d.test cyl3d.gpu.python.test cyl3d.mflow.test cyl3d.dtconst.test \ @@ -50,8 +51,8 @@ EXTRA_DIST = tap-driver.sh test_tps_splitcomm.py soln_differ inputs meshes lte- sgsSmag.test sgsSigma.test heatEq.test sponge.test plate.test pipe.mix.test lte2noneq-restart.test \ coupled-3d.interface.test plasma.axisym.test plasma.axisym.lte1d.test \ lomach-flow.test lomach-lequere.test interpInlet.test sgsLoMach.test autoPeriodic.test aveLoMach.test \ - reactFlow-binDiff.test reactFlow-singleRx.test reactFlow-table.test - + reactFlow-binDiff.test reactFlow-singleRx.test reactFlow-table.test \ + test_flow1d_coupling.py flow1d-coupling.python.test TESTS = vpath.sh XFAIL_TESTS = @@ -129,6 +130,11 @@ if PYTHON_ENABLED TESTS += cyl3d.python.test \ cyl3d.python.splitcomm.test \ coupled-3d.py-loop.test + +if GSLIB_ENABLED +TESTS += flow1d-coupling.python.test +endif + endif endif diff --git a/test/flow1d-coupling.python.test b/test/flow1d-coupling.python.test new file mode 100755 index 000000000..772251063 --- /dev/null +++ b/test/flow1d-coupling.python.test @@ -0,0 +1,41 @@ +#!./bats +# -*- mode: sh -*- + +TEST="flow1d-coupling" +RUNFILE="inputs/qms-axisym.flow1dcoupling.ini" +EXE="./test_flow1d_coupling.py" + +setup() { + SOLN_FILE=flow1d_coupling.sol.h5 + REF_FILE=ref_solns/flow1d_coupling.h5 +} + +@test "[$TEST] check for input file $RUNFILE" { + test -s $RUNFILE +} + +@test "[$TEST] run interface with input -> $RUNFILE" { + rm -f $SOLN_FILE + $EXE --runFile $RUNFILE + test -s $SOLN_FILE +} + +@test "[$TEST] verify interface output with input -> $RUNFILE" { + test -s $SOLN_FILE + test -s $REF_FILE + h5diff --report --relative=1e-12 $SOLN_FILE $REF_FILE /output/joule_heating +} + +@test "[$TEST] verify Joule heating is read only" { + # Should raise error + run WRITE_JOULE=True $EXE --runFile $RUNFILE + [ "$status" -ne 0 ] +} + +@test "[$TEST] verify total power input is 20 kW with input -> $RUNFILE { + TOTAL_POWER=True $EXE --runFile $RUNFILE +} + +@test "[$TEST] verify 1d and 2d total power inputs are equal with input -> $RUNFILE { + COMPARE_JOULE=True $EXE --runFile $RUNFILE +} \ No newline at end of file diff --git a/test/inputs/qms-axisym.flow1dcoupling.ini b/test/inputs/qms-axisym.flow1dcoupling.ini new file mode 100644 index 000000000..b653a627b --- /dev/null +++ b/test/inputs/qms-axisym.flow1dcoupling.ini @@ -0,0 +1,26 @@ +#--------------------- +# TPS runtime controls +#--------------------- + +# choice of solver +[solver] +type = em-axi # options are (flow, em, em-axi, or coupled) + +# controls for standalone EM-only solver +[em] +mesh = meshes/torch-em-coarse2.msh # mesh filename +order = 2 # FE order (polynomial degree) +max_iter = 60 # max number of iterations +rtol = 1.0e-13 # solver relative tolerance +atol = 1.0e-15 # solver absolute tolerance +top_only = false # run current through top branch only +bot_only = false # run current through bottom branch only +yinterp_min = -2.0 # minimum y interpolation value +yinterp_max = 2.0 # maximum y interpolation value +nBy = 129 # of interpolation points +By_file = ref_solns/By.h5 # file for By interpolant output + +# chosen for ~20 kW power input +current_amplitude = 4.575e7 # A/m^2 +current_frequency = 6e6 # 1/s +permeability = 1.25663706e-6 # m * kg / s^2 / A^2 \ No newline at end of file diff --git a/test/ref_solns/flow1d_coupling.h5 b/test/ref_solns/flow1d_coupling.h5 new file mode 100644 index 000000000..6c1dbd717 Binary files /dev/null and b/test/ref_solns/flow1d_coupling.h5 differ diff --git a/test/test_flow1d_coupling.py b/test/test_flow1d_coupling.py new file mode 100755 index 000000000..36ac916b8 --- /dev/null +++ b/test/test_flow1d_coupling.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 +""" This is a driver to test the interface between a Python 1d flow solver + and the 2d axisymmetric EM solver. The 1d flow solver is emulated by defining + the plasma conductivity""" +import sys +import os +import numpy as np +import h5py +from mpi4py import MPI +# set path to C++ TPS library +path = os.path.abspath(os.path.dirname(sys.argv[0])) +sys.path.append(path + "/../src/.libs") +import libtps + +comm = MPI.COMM_WORLD +# TPS solver +tps = libtps.Tps(comm) +tps.parseCommandLineArgs(sys.argv) +tps.parseInput() +tps.chooseDevices() + +# 1d flow - 2d EM interface +interface = libtps.Qms2Flow1d(tps) + +# Initialize MFEM vectors with number of 1d points +N_POINTS = 1000 +interface.initialize(N_POINTS) + +# Access vectors through NumPy +cond_1d = np.array(interface.PlasmaConductivity1d(), copy = False) +joule_1d = np.array(interface.JouleHeating1d(), copy = False) +radius_1d = np.array(interface.TorchRadius1d(), copy = False) +z_1d = np.array(interface.Coordinates1d(), copy = False) + +# Check that Joule heating is read only +if 'WRITE_JOULE' in os.environ: + joule_1d[0] = 1 + +# Example vector data +R_TORCH = 2.75e-2 +L_TORCH = 0.315 +z_1d[0:N_POINTS] = np.linspace(0, L_TORCH, N_POINTS) +cond_1d[0:N_POINTS] = np.ones(N_POINTS) +radius_1d[0:N_POINTS] = R_TORCH + +# Solve +N_INTERP = 100 +interface.set_n_interp(N_INTERP) +interface.solve() + +if 'TOTAL_POWER' in os.environ: + power = interface.total_joule_2d()/1e3 + assert np.abs(power - 20)/20 < 0.05, 'Total power should be 20 kW' + +if 'COMPARE_JOULE' in os.environ: + power_1d = interface.total_joule_1d()/1e3 + power_2d = interface.total_joule_2d()/1e3 + error = np.abs(power_1d - power_2d)/np.abs(power_2d) + assert error < 1e-5, '1d and 2d Joule heating integrals should match' + +SOLN_FILE = "flow1d_coupling.sol.h5" +# Save output with hdf5 +with h5py.File(SOLN_FILE, "w") as f: + _ = f.create_dataset('input/axial_coordinates', data=z_1d) + _ = f.create_dataset('input/torch_radius', data=radius_1d) + _ = f.create_dataset('input/plasma_conductivity', data=cond_1d) + _ = f.create_dataset('output/joule_heating', data=joule_1d) diff --git a/test/vpath.sh b/test/vpath.sh index af1e5eead..c4557e6e2 100755 --- a/test/vpath.sh +++ b/test/vpath.sh @@ -32,6 +32,7 @@ fi # necessary binaries binaries="bats die.sh soln_differ count_gpus.sh sniff_mpirun.sh " binaries+="../src/tps.py ../src/tps-time-loop.py ../src/tps-bte_0d3v.py ../test/test_tps_splitcomm.py" +binaries+=" ../test/test_flow1d_coupling.py" for binary in $binaries; do if [ ! -x $binary ];then if [ -x $testDir/$binary ];then