From 6e83ec087ab25cb1e3b94aace849635bea4cb1b7 Mon Sep 17 00:00:00 2001 From: Joshua Bigman Date: Thu, 25 Jul 2024 16:57:04 -0500 Subject: [PATCH 1/8] Created a new file for the CPUData functions which allow MFEM vectors to be read/write from Python as NumPy arrays. --- src/Makefile.am | 4 +-- src/data_exchange_utils.cpp | 68 +++++++++++++++++++++++++++++++++++++ src/data_exchange_utils.hpp | 67 ++++++++++++++++++++++++++++++++++++ src/pybindings.cpp | 4 +++ src/tps2Boltzmann.cpp | 49 +------------------------- 5 files changed, 142 insertions(+), 50 deletions(-) create mode 100644 src/data_exchange_utils.cpp create mode 100644 src/data_exchange_utils.hpp diff --git a/src/Makefile.am b/src/Makefile.am index 423a717d5..6cd546816 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -37,7 +37,7 @@ cxx_sources = averaging.cpp faceGradientIntegration.cpp M2ulPhyS.cpp rh 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 \ + 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,7 +53,7 @@ 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 \ + lte_transport_properties.hpp mixing_length_transport.hpp data_exchange_utils.hpp tps2Boltzmann.hpp gslib_interpolator.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 \ 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/pybindings.cpp b/src/pybindings.cpp index 2377b127c..c08cc2642 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/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) From 5f5a3c725885e6d08055bd33b8d0d9a8f62cef0c Mon Sep 17 00:00:00 2001 From: Joshua Bigman Date: Mon, 5 Aug 2024 11:19:38 -0500 Subject: [PATCH 2/8] Beginning of work on 1d flow-2d EM interface. The 2d plasma conductivity is set from the 1d values based on an assumed radial profile. --- src/Makefile.am | 4 +- src/flow1d-coupling.py | 33 +++++++ src/pybindings.cpp | 4 +- src/qms2flow1d.cpp | 200 +++++++++++++++++++++++++++++++++++++++++ src/qms2flow1d.hpp | 96 ++++++++++++++++++++ 5 files changed, 333 insertions(+), 4 deletions(-) create mode 100644 src/flow1d-coupling.py create mode 100644 src/qms2flow1d.cpp create mode 100644 src/qms2flow1d.hpp diff --git a/src/Makefile.am b/src/Makefile.am index 6cd546816..eee7cf4cb 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -36,7 +36,7 @@ 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 \ + 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 \ @@ -54,7 +54,7 @@ headers = averaging.hpp equation_of_state.hpp transport_properties. 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 data_exchange_utils.hpp tps2Boltzmann.hpp gslib_interpolator.hpp \ - loMach.hpp loMach_options.hpp calorically_perfect.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/flow1d-coupling.py b/src/flow1d-coupling.py new file mode 100644 index 000000000..47ac1fa1e --- /dev/null +++ b/src/flow1d-coupling.py @@ -0,0 +1,33 @@ +#!/usr/bin/env python3 +import sys +import os +import numpy as np +from mpi4py import MPI +# set path to C++ TPS library +path = os.path.abspath(os.path.dirname(sys.argv[0])) +sys.path.append(path + "/.libs") +print(sys.path) +import libtps + +comm = MPI.COMM_WORLD +# TPS solver +tps = libtps.Tps(comm) + +tps.parseCommandLineArgs(sys.argv) +tps.parseInput() +tps.chooseDevices() +#tps.chooseSolver() +#tps.initialize() + +print('Creating 1d interface instance') +interface = libtps.Qms2Flow1d(tps) +print('Initializing 1d interface vector size') +interface.initialize(3) +interface.print_all() +print('Writing plasma conductivity from Python') +cond_1d = np.array(interface.PlasmaConductivity1d(), copy = False) +cond_1d[0] = 1 +cond_1d[1] = 2 +cond_1d[2] = np.pi + +interface.print_all() diff --git a/src/pybindings.cpp b/src/pybindings.cpp index c08cc2642..074d0821d 100644 --- a/src/pybindings.cpp +++ b/src/pybindings.cpp @@ -44,7 +44,7 @@ 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); +void qms2flow1d(py::module& m); } // namespace tps_wrappers PYBIND11_MODULE(libtps, m) { @@ -53,7 +53,7 @@ PYBIND11_MODULE(libtps, m) { tps_wrappers::tps(m); tps_wrappers::tps2bolzmann(m); tps_wrappers::data_exchange_utils(m); - // tps_wrappers::qms2flow1d(m); + tps_wrappers::qms2flow1d(m); } #endif diff --git a/src/qms2flow1d.cpp b/src/qms2flow1d.cpp new file mode 100644 index 000000000..f7680d4b6 --- /dev/null +++ b/src/qms2flow1d.cpp @@ -0,0 +1,200 @@ +// -----------------------------------------------------------------------------------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; + + n_1d = 0; + cond_1d = nullptr; + joule_1d = nullptr; + z_coords_1d = nullptr; + radius_1d = nullptr; +} + +Qms2Flow1d::~Qms2Flow1d() { + delete cond_1d; + delete joule_1d; + delete z_coords_1d; + delete radius_1d; + + delete em_opts; + delete qmsa; +} + +void Qms2Flow1d::initialize(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_); + + qmsa->parseSolverOptions(); + qmsa->initialize(); +} + +void Qms2Flow1d::print_all_1d() { + std::cout << "---\n1d Coordinates:\n---" << endl; + z_coords_1d->Print(); + std::cout << "---\nTorch Radius:\n---" << endl; + radius_1d->Print(); + std::cout << "---\n1d Plasma Conductivity:\n---" << endl; + cond_1d->Print(); + std::cout << "---\n1d Joule Heating:\n---" << endl; + joule_1d->Print(); +} + +double Qms2Flow1d::radial_profile(double r, double R, double r_c) const { + // Constant over [0, r_c]. Zero beyond torch radius R + // Cosine for cont. connection between the two + assert(0 <= r && 0 <= R && 0 <= r_c); + double out = 0; + if (0 <= r && r <= r_c) { + out = 1; + } else if (r_c < r && r < R) { + out = 1 + cos(M_PI*(r - r_c)/(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.9*torch_radius); + sigma *= sigma_1d; + return sigma; +} + +void Qms2Flow1d::set_plasma_conductivity() { + ParFiniteElementSpace *fes = qmsa->getFESpace(); + ParMesh *pmesh = qmsa->getMesh(); + + ParGridFunction *coordsDof = new ParGridFunction(fes); + pmesh->GetNodes(*coordsDof); + + ParGridFunction *pc = qmsa->getPlasmaConductivityGF(); + double *plasma_conductivity_gf = pc->HostWrite(); + + const int n_nodes = pc->FESpace()->GetNDofs(); + for (int i = 0; i < n_nodes; i++) { + const double r = (*coordsDof)[i]; + const double z = (*coordsDof)[i + n_nodes]; + plasma_conductivity_gf[i] = expand_cond_2d(r, z); + } + + delete coordsDof; +} + +} // 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) + .def("print_all", &TPS::Qms2Flow1d::print_all_1d) + .def("Coordinates1d", + [](TPS::Qms2Flow1d &interface) { + return std::unique_ptr(new TPS::CPUData(interface.Coordinates1d(), false)); + }) + .def("Radius1d", + [](TPS::Qms2Flow1d &interface) { + return std::unique_ptr(new TPS::CPUData(interface.Radius1d(), false)); + }) + .def("PlasmaConductivity1d", + [](TPS::Qms2Flow1d &interface) { + return std::unique_ptr(new TPS::CPUData(interface.PlasmaConductivity1d(), false)); + }) + .def("JouleHeating1d", + [](TPS::Qms2Flow1d &interface) { + return std::unique_ptr(new TPS::CPUDataRead(interface.JouleHeating1d())); + }); +} + +} // namespace tps_wrappers + +#endif diff --git a/src/qms2flow1d.hpp b/src/qms2flow1d.hpp new file mode 100644 index 000000000..8d5993942 --- /dev/null +++ b/src/qms2flow1d.hpp @@ -0,0 +1,96 @@ +// -----------------------------------------------------------------------------------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" + +namespace TPS { + +// Need to: +// [X] - Initialize vector of integer size with argument given in Python +// (Number of 1D grid points) +// [X] - Form vector of that size in C++ +// (To receive plasma conductivity/send Joule heating) +// [X] - Receive 1D plasma conductivity from Python +// [ ] - Expand to 2D plasma conductivity in C++ +// [ ] - Solve for 2D Joule heating with axisymmetric solver +// [ ] - Reduce to 1D Joule heating in C++ +// [X] - Send 1D Joule heating to Python + +class Qms2Flow1d { + public: + Qms2Flow1d(Tps *tps); + ~Qms2Flow1d(); + void initialize(int n_1d_); + 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 &Radius1d() const { return *radius_1d; } + Vector &Radius1d() { return *radius_1d; } + + private: + Tps *tpsP_; + ElectromagneticOptions *em_opts; + QuasiMagnetostaticSolverAxiSym *qmsa; + + int n_1d; + Vector *cond_1d; + Vector *joule_1d; + Vector *z_coords_1d; // 1d coordinates are assumed to be sorted + Vector *radius_1d; // Torch radius at z coordinates + + // r is radial location, R is torch radius, r_c is modeling parameter + double radial_profile(double r, double R, double r_c) const; + // Binary search to find interval + int find_z_interval(double z) const; + // Interpolation of 1d values along centerline + double interpolate_z(double z, const Vector &values) const; + double expand_cond_2d(double r, double z) const; + + void set_plasma_conductivity(); +}; + +} // namespace TPS +#endif // QMS2FLOW1D_HPP_ From d8de08d5a936457c3fd0c5307184a4cdd8e1800b Mon Sep 17 00:00:00 2001 From: Joshua Bigman Date: Wed, 7 Aug 2024 21:53:10 -0500 Subject: [PATCH 3/8] Created LineInterpolator class using GSLib. Fixed a segfault issue when setting the 2d plasma conductivity and created a Python function to update the plasma conductivity then solve the EM field. --- src/flow1d-coupling.py | 25 ++++++++++++++------- src/gslib_interpolator.cpp | 45 ++++++++++++++++++++++++++++++++++++++ src/gslib_interpolator.hpp | 22 +++++++++++++++++++ src/qms2flow1d.cpp | 21 +++++++++++++----- src/qms2flow1d.hpp | 8 ++++--- 5 files changed, 104 insertions(+), 17 deletions(-) diff --git a/src/flow1d-coupling.py b/src/flow1d-coupling.py index 47ac1fa1e..b78ca9c92 100644 --- a/src/flow1d-coupling.py +++ b/src/flow1d-coupling.py @@ -6,7 +6,6 @@ # set path to C++ TPS library path = os.path.abspath(os.path.dirname(sys.argv[0])) sys.path.append(path + "/.libs") -print(sys.path) import libtps comm = MPI.COMM_WORLD @@ -19,15 +18,25 @@ #tps.chooseSolver() #tps.initialize() -print('Creating 1d interface instance') +print('-'*40, '\nCreating 1d interface instance\n', '-'*40) interface = libtps.Qms2Flow1d(tps) -print('Initializing 1d interface vector size') -interface.initialize(3) +print('-'*40, '\nInitializing 1d interface vector size\n', '-'*40) +n_points = 25 +interface.initialize(n_points) interface.print_all() -print('Writing plasma conductivity from Python') + +print('-'*40, '\nCreating NumPy array interfaces for vectors\n', '-'*40) cond_1d = np.array(interface.PlasmaConductivity1d(), copy = False) -cond_1d[0] = 1 -cond_1d[1] = 2 -cond_1d[2] = np.pi +joule_1d = np.array(interface.JouleHeating1d(), copy = False) +radius_1d = np.array(interface.TorchRadius1d(), copy = False) +z_1d = np.array(interface.Coordinates1d(), copy = False) +R_TORCH = 2.75e-2 +print('-'*40, '\nSetting vector values\n', '-'*40) +z_1d[0:n_points] = np.linspace(0, 0.5, n_points) +cond_1d[0:n_points] = np.ones(n_points) +radius_1d[0:n_points] = R_TORCH*(1 - 0.1*np.exp(-((z_1d - 0.25)/0.05)**2)) interface.print_all() + +print('-'*40, '\nSolving EM field\n', '-'*40) +interface.solve() diff --git a/src/gslib_interpolator.cpp b/src/gslib_interpolator.cpp index 5046e57e6..e95e379d2 100644 --- a/src/gslib_interpolator.cpp +++ b/src/gslib_interpolator.cpp @@ -204,3 +204,48 @@ 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() : InterpolatorBase() { + 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::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..fcbb7cc65 100644 --- a/src/gslib_interpolator.hpp +++ b/src/gslib_interpolator.hpp @@ -110,4 +110,26 @@ 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(); + LineInterpolator(mfem::Vector start, mfem::Vector end, int n); + virtual ~LineInterpolator(); + + /** 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/qms2flow1d.cpp b/src/qms2flow1d.cpp index f7680d4b6..4e19aa89c 100644 --- a/src/qms2flow1d.cpp +++ b/src/qms2flow1d.cpp @@ -81,6 +81,11 @@ void Qms2Flow1d::initialize(int n_1d_) { qmsa->initialize(); } +void Qms2Flow1d::solve() { + set_plasma_conductivity(); + qmsa->solve(); +} + void Qms2Flow1d::print_all_1d() { std::cout << "---\n1d Coordinates:\n---" << endl; z_coords_1d->Print(); @@ -150,19 +155,22 @@ void Qms2Flow1d::set_plasma_conductivity() { ParFiniteElementSpace *fes = qmsa->getFESpace(); ParMesh *pmesh = qmsa->getMesh(); - ParGridFunction *coordsDof = new ParGridFunction(fes); + 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 = pc->FESpace()->GetNDofs(); + const int n_nodes = fes->GetNDofs(); + double r, z; for (int i = 0; i < n_nodes; i++) { - const double r = (*coordsDof)[i]; - const double z = (*coordsDof)[i + n_nodes]; + r = coordsDof->Elem(i); + z = coordsDof->Elem(i + n_nodes); plasma_conductivity_gf[i] = expand_cond_2d(r, z); } + delete coord_fes; delete coordsDof; } @@ -176,14 +184,15 @@ void qms2flow1d(py::module &m) { py::class_(m, "Qms2Flow1d") .def(py::init()) .def("initialize", &TPS::Qms2Flow1d::initialize) + .def("solve", &TPS::Qms2Flow1d::solve) .def("print_all", &TPS::Qms2Flow1d::print_all_1d) .def("Coordinates1d", [](TPS::Qms2Flow1d &interface) { return std::unique_ptr(new TPS::CPUData(interface.Coordinates1d(), false)); }) - .def("Radius1d", + .def("TorchRadius1d", [](TPS::Qms2Flow1d &interface) { - return std::unique_ptr(new TPS::CPUData(interface.Radius1d(), false)); + return std::unique_ptr(new TPS::CPUData(interface.TorchRadius1d(), false)); }) .def("PlasmaConductivity1d", [](TPS::Qms2Flow1d &interface) { diff --git a/src/qms2flow1d.hpp b/src/qms2flow1d.hpp index 8d5993942..f1e117ef7 100644 --- a/src/qms2flow1d.hpp +++ b/src/qms2flow1d.hpp @@ -46,7 +46,7 @@ namespace TPS { // [X] - Form vector of that size in C++ // (To receive plasma conductivity/send Joule heating) // [X] - Receive 1D plasma conductivity from Python -// [ ] - Expand to 2D plasma conductivity in C++ +// [X] - Expand to 2D plasma conductivity in C++ // [ ] - Solve for 2D Joule heating with axisymmetric solver // [ ] - Reduce to 1D Joule heating in C++ // [X] - Send 1D Joule heating to Python @@ -56,6 +56,7 @@ class Qms2Flow1d { Qms2Flow1d(Tps *tps); ~Qms2Flow1d(); void initialize(int n_1d_); + void solve(); // Set the plasma conductivity and solve the EM field void print_all_1d(); const Vector &PlasmaConductivity1d() const { return *cond_1d; } @@ -67,8 +68,8 @@ class Qms2Flow1d { const Vector &Coordinates1d() const { return *z_coords_1d; } Vector &Coordinates1d() { return *z_coords_1d; } - const Vector &Radius1d() const { return *radius_1d; } - Vector &Radius1d() { return *radius_1d; } + const Vector &TorchRadius1d() const { return *radius_1d; } + Vector &TorchRadius1d() { return *radius_1d; } private: Tps *tpsP_; @@ -87,6 +88,7 @@ class Qms2Flow1d { int find_z_interval(double z) const; // Interpolation of 1d values along centerline double interpolate_z(double z, const Vector &values) const; + double expand_cond_2d(double r, double z) const; void set_plasma_conductivity(); From 4e832ee60ebcfb6b39d2fe9ad59b93bbb73a677e Mon Sep 17 00:00:00 2001 From: Joshua Bigman Date: Fri, 9 Aug 2024 17:00:13 -0500 Subject: [PATCH 4/8] Added 1d Joule heating calculations based on interpolation and radial integration. --- src/flow1d-coupling.py | 6 ++-- src/gslib_interpolator.cpp | 11 ++++++ src/gslib_interpolator.hpp | 10 ++++++ src/qms2flow1d.cpp | 71 +++++++++++++++++++++++++++++++++----- src/qms2flow1d.hpp | 24 ++++++------- 5 files changed, 96 insertions(+), 26 deletions(-) diff --git a/src/flow1d-coupling.py b/src/flow1d-coupling.py index b78ca9c92..f22780f3c 100644 --- a/src/flow1d-coupling.py +++ b/src/flow1d-coupling.py @@ -11,12 +11,9 @@ comm = MPI.COMM_WORLD # TPS solver tps = libtps.Tps(comm) - tps.parseCommandLineArgs(sys.argv) tps.parseInput() tps.chooseDevices() -#tps.chooseSolver() -#tps.initialize() print('-'*40, '\nCreating 1d interface instance\n', '-'*40) interface = libtps.Qms2Flow1d(tps) @@ -31,7 +28,7 @@ radius_1d = np.array(interface.TorchRadius1d(), copy = False) z_1d = np.array(interface.Coordinates1d(), copy = False) -R_TORCH = 2.75e-2 +R_TORCH = 2.75e-1 print('-'*40, '\nSetting vector values\n', '-'*40) z_1d[0:n_points] = np.linspace(0, 0.5, n_points) cond_1d[0:n_points] = np.ones(n_points) @@ -40,3 +37,4 @@ print('-'*40, '\nSolving EM field\n', '-'*40) interface.solve() +interface.print_all() diff --git a/src/gslib_interpolator.cpp b/src/gslib_interpolator.cpp index e95e379d2..5f8f3e531 100644 --- a/src/gslib_interpolator.cpp +++ b/src/gslib_interpolator.cpp @@ -210,11 +210,22 @@ LineInterpolator::LineInterpolator() : InterpolatorBase() { end_.SetSize(0); } +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_); diff --git a/src/gslib_interpolator.hpp b/src/gslib_interpolator.hpp index fcbb7cc65..f482e77cf 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 { @@ -121,14 +124,21 @@ class LineInterpolator : public InterpolatorBase { public: LineInterpolator(); + LineInterpolator(int n); LineInterpolator(mfem::Vector start, mfem::Vector end, int n); virtual ~LineInterpolator(); + void initializeFinder(mfem::ParMesh *mesh); + + void setStart(mfem::Vector start) { start_ = start; } + void setEnd(mfem::Vector end) { end_ = 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; }; diff --git a/src/qms2flow1d.cpp b/src/qms2flow1d.cpp index 4e19aa89c..965493c84 100644 --- a/src/qms2flow1d.cpp +++ b/src/qms2flow1d.cpp @@ -32,6 +32,7 @@ #include "qms2flow1d.hpp" #include "data_exchange_utils.hpp" +#include "gslib_interpolator.hpp" #ifdef HAVE_PYTHON #include @@ -79,11 +80,14 @@ void Qms2Flow1d::initialize(int n_1d_) { qmsa->parseSolverOptions(); qmsa->initialize(); + // Ensuring nodes is required for GSLib interpolation + qmsa->getMesh()->EnsureNodes(); } void Qms2Flow1d::solve() { - set_plasma_conductivity(); + set_plasma_conductivity_2d(); qmsa->solve(); + set_joule_heating_1d(10); } void Qms2Flow1d::print_all_1d() { @@ -97,15 +101,15 @@ void Qms2Flow1d::print_all_1d() { joule_1d->Print(); } -double Qms2Flow1d::radial_profile(double r, double R, double r_c) const { - // Constant over [0, r_c]. Zero beyond torch radius R +double Qms2Flow1d::radial_profile(double r, double torch_r, double r_c) const { + // Constant over [0, r_c]. Zero beyond torch_r // Cosine for cont. connection between the two - assert(0 <= r && 0 <= R && 0 <= r_c); + 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 < R) { - out = 1 + cos(M_PI*(r - r_c)/(R - r_c)); + } 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; @@ -146,12 +150,12 @@ double Qms2Flow1d::interpolate_z(double z, const Vector &values) const { 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.9*torch_radius); + double sigma = radial_profile(r, torch_radius, 0.7*torch_radius); sigma *= sigma_1d; return sigma; } -void Qms2Flow1d::set_plasma_conductivity() { +void Qms2Flow1d::set_plasma_conductivity_2d() { ParFiniteElementSpace *fes = qmsa->getFESpace(); ParMesh *pmesh = qmsa->getMesh(); @@ -174,6 +178,57 @@ void Qms2Flow1d::set_plasma_conductivity() { delete coordsDof; } +// Trapezoidal rule +double radial_integral(Vector &values, Vector &r_coords) { + assert(values.Size() == r_coords.Size()); + const int n_minus_one = values.Size() - 1; + + // r_i - r_i-1 + Vector diff = Vector(r_coords, 0, n_minus_one); + diff -= Vector(r_coords, 1, n_minus_one); + + // v_i*r_i + values *= r_coords; + + // (v_i*r_i + v_i-1*r_i-1)(r_i - r_i-1) + Vector sum = Vector(values, 0, n_minus_one); + sum += Vector(values, 1, n_minus_one); + sum *= diff; + + // Factor of 2pi for angles and division by 2 in trap. rule cancel + double integral = M_PI*sum.Sum(); + + return integral; +} + +void Qms2Flow1d::set_joule_heating_1d(const int n_interp) { + const int dim = 2; + double z, torch_r; + Vector start = Vector(dim), end = Vector(dim); + LineInterpolator l_interp = LineInterpolator(n_interp); + l_interp.initializeFinder(qmsa->getMesh()); + + for (int i = 0; i < n_1d; i++) { + z = z_coords_1d->Elem(i); + torch_r = radius_1d->Elem(i); + + start[0] = 0; + start[1] = z; + end[0] = torch_r; + end[1] = z; + + l_interp.setStart(start); + l_interp.setEnd(end); + 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 diff --git a/src/qms2flow1d.hpp b/src/qms2flow1d.hpp index f1e117ef7..85a2e1853 100644 --- a/src/qms2flow1d.hpp +++ b/src/qms2flow1d.hpp @@ -40,17 +40,6 @@ namespace TPS { -// Need to: -// [X] - Initialize vector of integer size with argument given in Python -// (Number of 1D grid points) -// [X] - Form vector of that size in C++ -// (To receive plasma conductivity/send Joule heating) -// [X] - Receive 1D plasma conductivity from Python -// [X] - Expand to 2D plasma conductivity in C++ -// [ ] - Solve for 2D Joule heating with axisymmetric solver -// [ ] - Reduce to 1D Joule heating in C++ -// [X] - Send 1D Joule heating to Python - class Qms2Flow1d { public: Qms2Flow1d(Tps *tps); @@ -82,8 +71,12 @@ class Qms2Flow1d { Vector *z_coords_1d; // 1d coordinates are assumed to be sorted Vector *radius_1d; // Torch radius at z coordinates - // r is radial location, R is torch radius, r_c is modeling parameter - double radial_profile(double r, double R, double r_c) const; + /** + * 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 @@ -91,7 +84,10 @@ class Qms2Flow1d { double expand_cond_2d(double r, double z) const; - void set_plasma_conductivity(); + void set_plasma_conductivity_2d(); + + int n_interp; + void set_joule_heating_1d(const int n_interp); }; } // namespace TPS From 5177edb691b14fbd219a2f10e1778968a2f67ebc Mon Sep 17 00:00:00 2001 From: Joshua Bigman Date: Mon, 12 Aug 2024 15:44:28 -0500 Subject: [PATCH 5/8] Separated numerical integrals into linear and radial, and fixed an overwriting error in the implementation. Changed how the LineInterpolator is implemented in the interface and fixed a memory leak due to internal data remaining in FindPointsGSLib. Added additional methods to control the interface behavior from Python. Added various comments and descriptions to methods visible from Python. --- src/flow1d-coupling.py | 31 +++++---- src/gslib_interpolator.cpp | 7 +- src/gslib_interpolator.hpp | 6 +- src/qms2flow1d.cpp | 133 ++++++++++++++++++++++--------------- src/qms2flow1d.hpp | 55 ++++++++++----- 5 files changed, 143 insertions(+), 89 deletions(-) diff --git a/src/flow1d-coupling.py b/src/flow1d-coupling.py index f22780f3c..c4801c937 100644 --- a/src/flow1d-coupling.py +++ b/src/flow1d-coupling.py @@ -15,26 +15,33 @@ tps.parseInput() tps.chooseDevices() -print('-'*40, '\nCreating 1d interface instance\n', '-'*40) +# 1d flow - 2d EM interface interface = libtps.Qms2Flow1d(tps) -print('-'*40, '\nInitializing 1d interface vector size\n', '-'*40) -n_points = 25 +# Initialize MFEM vectors with number of 1d points +n_points = 1000 interface.initialize(n_points) -interface.print_all() -print('-'*40, '\nCreating NumPy array interfaces for vectors\n', '-'*40) +# 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) -R_TORCH = 2.75e-1 -print('-'*40, '\nSetting vector values\n', '-'*40) -z_1d[0:n_points] = np.linspace(0, 0.5, n_points) +# Example vector data +R_TORCH = 2.75e-2 +z_1d[0:n_points] = np.linspace(0, 0.315, n_points) cond_1d[0:n_points] = np.ones(n_points) -radius_1d[0:n_points] = R_TORCH*(1 - 0.1*np.exp(-((z_1d - 0.25)/0.05)**2)) -interface.print_all() +radius_1d[0:n_points] = R_TORCH -print('-'*40, '\nSolving EM field\n', '-'*40) +# Solve +interface.set_n_interp(100) interface.solve() -interface.print_all() + +# Evaluate Joule heating +tot_1d = interface.total_joule_1d()/1e3 +tot_2d = interface.total_joule_2d()/1e3 +error = 100*np.abs(tot_1d - tot_2d)/np.abs(tot_2d) + +print(f'Total 2d Joule heating [kW]: {tot_2d:.2f}') +print(f'Total 1d Joule heating [kW]: {tot_1d:.2f}') +print(f'Error [%]: {error:.3e}') diff --git a/src/gslib_interpolator.cpp b/src/gslib_interpolator.cpp index 5f8f3e531..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 } @@ -205,11 +207,6 @@ void PlaneInterpolator::writeAscii(std::string oname, bool rank0) const { InterpolatorBase::writeAscii(oname, rank0); } -LineInterpolator::LineInterpolator() : InterpolatorBase() { - start_.SetSize(0); - end_.SetSize(0); -} - LineInterpolator::LineInterpolator(int n) : InterpolatorBase(), n_(n) { start_.SetSize(0); end_.SetSize(0); diff --git a/src/gslib_interpolator.hpp b/src/gslib_interpolator.hpp index f482e77cf..d7d81364b 100644 --- a/src/gslib_interpolator.hpp +++ b/src/gslib_interpolator.hpp @@ -123,15 +123,17 @@ class LineInterpolator : public InterpolatorBase { int n_; public: - LineInterpolator(); - LineInterpolator(int n); + 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 * diff --git a/src/qms2flow1d.cpp b/src/qms2flow1d.cpp index 965493c84..f1d36e07e 100644 --- a/src/qms2flow1d.cpp +++ b/src/qms2flow1d.cpp @@ -32,7 +32,6 @@ #include "qms2flow1d.hpp" #include "data_exchange_utils.hpp" -#include "gslib_interpolator.hpp" #ifdef HAVE_PYTHON #include @@ -50,12 +49,14 @@ 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() { @@ -66,17 +67,21 @@ Qms2Flow1d::~Qms2Flow1d() { delete em_opts; delete qmsa; + delete l_interp; } -void Qms2Flow1d::initialize(int n_1d_) { +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_); + 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(); @@ -87,23 +92,23 @@ void Qms2Flow1d::initialize(int n_1d_) { void Qms2Flow1d::solve() { set_plasma_conductivity_2d(); qmsa->solve(); - set_joule_heating_1d(10); + set_joule_heating_1d(); } void Qms2Flow1d::print_all_1d() { - std::cout << "---\n1d Coordinates:\n---" << endl; + std::cout << "\n--- 1d Coordinates ---" << endl; z_coords_1d->Print(); - std::cout << "---\nTorch Radius:\n---" << endl; + std::cout << "\n--- Torch Radius ---" << endl; radius_1d->Print(); - std::cout << "---\n1d Plasma Conductivity:\n---" << endl; + std::cout << "\n--- 1d Plasma Conductivity ---" << endl; cond_1d->Print(); - std::cout << "---\n1d Joule Heating:\n---" << endl; + 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 cont. connection between the two + // 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) { @@ -178,54 +183,60 @@ void Qms2Flow1d::set_plasma_conductivity_2d() { delete coordsDof; } -// Trapezoidal rule -double radial_integral(Vector &values, Vector &r_coords) { - assert(values.Size() == r_coords.Size()); +double Qms2Flow1d::linear_integral(Vector values, Vector coords) { + assert(values.Size() == coords.Size()); const int n_minus_one = values.Size() - 1; - // r_i - r_i-1 - Vector diff = Vector(r_coords, 0, n_minus_one); - diff -= Vector(r_coords, 1, n_minus_one); - - // v_i*r_i - values *= r_coords; - - // (v_i*r_i + v_i-1*r_i-1)(r_i - r_i-1) - Vector sum = Vector(values, 0, n_minus_one); - sum += Vector(values, 1, n_minus_one); - sum *= diff; - - // Factor of 2pi for angles and division by 2 in trap. rule cancel - double integral = M_PI*sum.Sum(); + // 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(); +} - return integral; +// 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(const int n_interp) { - const int dim = 2; +void Qms2Flow1d::set_joule_heating_1d() { double z, torch_r; - Vector start = Vector(dim), end = Vector(dim); - LineInterpolator l_interp = LineInterpolator(n_interp); - l_interp.initializeFinder(qmsa->getMesh()); + l_interp->setn(n_interp); + l_interp->initializeFinder(qmsa->getMesh()); for (int i = 0; i < n_1d; i++) { z = z_coords_1d->Elem(i); - torch_r = radius_1d->Elem(i); - start[0] = 0; - start[1] = z; - end[0] = torch_r; - end[1] = z; + 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.setStart(start); - l_interp.setEnd(end); - l_interp.setInterpolationPoints(); - l_interp.interpolate(qmsa->getJouleHeatingGF()); + l_interp->setInterpolationPoints(); + l_interp->interpolate(qmsa->getJouleHeatingGF()); - Vector values = l_interp.getSolution(); - Vector r_coords = Vector(l_interp.getInterpolationPoints(), 0, n_interp); + Vector values = l_interp->getSolution(); + Vector r_coords = Vector(l_interp->getInterpolationPoints(), 0, n_interp); - joule_1d->Elem(i) = radial_integral(values, r_coords); + joule_1d->Elem(i) = radial_integral(values, r_coords); + } } } @@ -238,25 +249,39 @@ namespace tps_wrappers { void qms2flow1d(py::module &m) { py::class_(m, "Qms2Flow1d") .def(py::init()) - .def("initialize", &TPS::Qms2Flow1d::initialize) - .def("solve", &TPS::Qms2Flow1d::solve) - .def("print_all", &TPS::Qms2Flow1d::print_all_1d) + .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 diff --git a/src/qms2flow1d.hpp b/src/qms2flow1d.hpp index 85a2e1853..489878fc2 100644 --- a/src/qms2flow1d.hpp +++ b/src/qms2flow1d.hpp @@ -37,6 +37,7 @@ #include "tps.hpp" #include "em_options.hpp" #include "quasimagnetostatic.hpp" +#include "gslib_interpolator.hpp" namespace TPS { @@ -44,8 +45,11 @@ class Qms2Flow1d { public: Qms2Flow1d(Tps *tps); ~Qms2Flow1d(); - void initialize(int n_1d_); - void solve(); // Set the plasma conductivity and solve the EM field + // 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; } @@ -60,34 +64,53 @@ class Qms2Flow1d { 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; - Vector *z_coords_1d; // 1d coordinates are assumed to be sorted - Vector *radius_1d; // Torch radius at z coordinates - - /** - * r: radial location - * torch_r: torch radius - * r_c: modeling parameter, 0 <= r_c <= R - */ + // 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 + // 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(); - - int n_interp; - void set_joule_heating_1d(const int n_interp); }; } // namespace TPS From 27abd90d19c95cfffc72e29449299851a6c2ba14 Mon Sep 17 00:00:00 2001 From: Joshua Bigman Date: Thu, 22 Aug 2024 10:43:54 -0500 Subject: [PATCH 6/8] Created tests for the 1d flow - axisymmetric EM coupling. Tests verify the code executes without error and matches a reference solution. It does not check for accuracy to a known or manufactured solution. --- .gitattributes | 1 + docker/test/Dockerfile | 1 + src/flow1d-coupling.py | 47 --------------- test/Makefile.am | 10 +++- test/flow1d-coupling.python.test | 41 ++++++++++++++ test/inputs/qms-axisym.flow1dcoupling.ini | 26 +++++++++ test/ref_solns/flow1d_coupling.h5 | Bin 0 -> 36480 bytes test/test_flow1d_coupling.py | 66 ++++++++++++++++++++++ test/vpath.sh | 1 + 9 files changed, 144 insertions(+), 49 deletions(-) delete mode 100644 src/flow1d-coupling.py create mode 100755 test/flow1d-coupling.python.test create mode 100644 test/inputs/qms-axisym.flow1dcoupling.ini create mode 100644 test/ref_solns/flow1d_coupling.h5 create mode 100755 test/test_flow1d_coupling.py diff --git a/.gitattributes b/.gitattributes index dda8bcaaf..73e0aca5d 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.sol.h5 filter=lfs diff=lfs merge=lfs -text diff --git a/docker/test/Dockerfile b/docker/test/Dockerfile index d3869b629..3ec07450a 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 diff --git a/src/flow1d-coupling.py b/src/flow1d-coupling.py deleted file mode 100644 index c4801c937..000000000 --- a/src/flow1d-coupling.py +++ /dev/null @@ -1,47 +0,0 @@ -#!/usr/bin/env python3 -import sys -import os -import numpy as np -from mpi4py import MPI -# set path to C++ TPS library -path = os.path.abspath(os.path.dirname(sys.argv[0])) -sys.path.append(path + "/.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) - -# Example vector data -R_TORCH = 2.75e-2 -z_1d[0:n_points] = np.linspace(0, 0.315, n_points) -cond_1d[0:n_points] = np.ones(n_points) -radius_1d[0:n_points] = R_TORCH - -# Solve -interface.set_n_interp(100) -interface.solve() - -# Evaluate Joule heating -tot_1d = interface.total_joule_1d()/1e3 -tot_2d = interface.total_joule_2d()/1e3 -error = 100*np.abs(tot_1d - tot_2d)/np.abs(tot_2d) - -print(f'Total 2d Joule heating [kW]: {tot_2d:.2f}') -print(f'Total 1d Joule heating [kW]: {tot_1d:.2f}') -print(f'Error [%]: {error:.3e}') diff --git a/test/Makefile.am b/test/Makefile.am index 0366d3e99..5d0ef1593 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.sol.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 0000000000000000000000000000000000000000..6c1dbd7176cbf7f1947912bdd1142c36c99a40a1 GIT binary patch literal 36480 zcmeI5c{EkwzyHni(4a|!LWM$ur1AM|B@&quN=ivYWek}smwE(=Xstv15H#S zW0DLRqV73*x@+C<`rYqc_n-T_>w9UfwsZDA`|Rg=KF|C8dOv%wy03Kj2rri)7X$r= zgM)#EVbT1rfB*W8f9^i|^UmMbfA^omzcAvh2j0dp{rQEFf&PbqVI_WlF+Tr49(O`h zL4kpGuJaGN| zZKa{TsiFDV3)ik$8JS)+w7F=Fci3FBx?pnF%FxKvcK$63b3^MZhO`5(8rfd3F}-eT z<1l~Z-+KS6lUm2s4l^<^G2rrJSdZ)W?{WWF$A3T9De>ov{~pJ&<sOkNfrSmYKGnzGv}zg!>*CWVjh*e&gd~) zK0lpYlOTisY!4&dlP1<%ejmrLPax&zPP?L7_3--!DAPSw_pQKryx$0w!=r$6=n_88 z6p^QGTKi1R@$pxY#aqeyyX!6R=d94=t_@;0&syWp+n_tX0yV|sw)i}Di0vD@F_Wh~ zKA!{fI`)(0h|o=Z-dku{+;HUmOh)U>v2b)R~);Y_Sqd@mj@E7 zu2j{MzK^f(J~FpCxOT$96JO^8lp>rv6j1p9U++USM%rE|$M1!&+Y9Z!llb7B-Xnbd zkI)Utji+x#dE<6GMzo7r=56;Fx92fpRN$HWxW)&!%Lj?9sr$UG+y}SM7b%RnH9B4L z#qIP(#`W@<6rLxzy-$!w`0o&r^e4F8ekj4d(8Wf_54YbBHL5w7Rn7Y0e(*=bn`E_? zNBHCZ@JEY&jqbgq7=ZgF0BvuZZz1lVK%|vq&c*2)i2ErJ*}902>g^7~{S|})PxZb{ z=m^677L4+?N0$9`3da2xjM_QQ#j9+7iu>^?qBr(=<5S$9A!uFBTl-y5yUVu9u3! z;~s^qUnUmS-i^ZJABC6&J-1s@qi`KWAwM^h=vAAdaXmz%jo*)~Dz}fubrFrS^*4H5 zZivS95segUnHiU^i@|jggWl}yvMI2L!Sxb@OhZfiFI2?fx`{zET%r0r%VTl<#G(h+ z9Z8ucvAB+6(VDK)htCwm;(Cfj$*OyES^45{UB#h;c`LWRJQs)SD-O{I>tse8uCq9F z_HomVt-r~{7&kUMk#pAk*N3^-}Q!Ozb*IztZLbbRyH5ZTTFagaEuX+Nm#{?vk zAgOaSG6C0R0;*cZ)9`&F0oP|DI_Wfc@TpQFuG2)+H&mN0ACidcH4)kC#%}uZH4)cs z65_0Md-?>DaQ!Bsr_x4z`=2D?I!;1cLgbx4_axzZPDTZ68-6^LOU8Abj8v?Co9=p$ zjO#lYy`TR{#&!M@T|HW;ac|p8T<VS*E{Ta-S!gK|4X#tp-%E0 zu@pQHQc(K5K&RmOkb(|f5b!58rQmsyf-YZjpiQw#=j^`%DkN2eFd6J5LQgo(1 zpHIc}B^6b0DDK$vFcr_6RP^wmthrZKDxN>7=!DJ2m`@*4@jOaJYZC-IcX6cQ`ILrw zKGG{R4bQ7I^m5sl)`wGRcz&g!>&LphcDSbDd6tF_I+1JcC8y!}mWJ4KN*UW3CkIqwf{wB6p+H z@%&6j`O~M*v^J;Xd76%#rH=)QO{L@cnt_g8+)KW-Ap_6b474;voSUqif#+`q>L9ID z*l3f1=Wzy#V&k~u@H_+0=L}@Ne=@hJIs?z^47AItcXY!@2A-ANcc%EmX z>i3m&3&8U|6Fpj*8Me+m6VLlhL@VZZR{oiI{%0bgc>~VGI*^I_Go6mF?$5+}kcsB) zba58eg)HQ-{#<#*-Yl#SS%^5PHoL+o3+qG{;`%edSTC{=rA*SIG&>9HMi!#ehCpW) z){iW7L6YaQ31>Fek!-Z(!q2ruJF>B!WTWw*L7m0=*;rSyQDO62pNnqUSYNV{8%te7 z-pg#PGudcfF5YBgy~#$)En?Am#vH6WIp{-lfJ3%K4%VL>6x;5W!F3`B>rf85w%GQo zfnyHVqZ}lsXtX&aCI{5o{M!Y4;|jpv@KpR z59?bV;xR0|Ja;G$>s%gc3rLL8xRQtUE)UVXl-A@ttb2Kg23IR*O7gJ&yO@vlF(0Y(%kB+*l#g{XAJMCM zEGHl9Wj^|JL*S+2r+loN`6xb>xj%%f0PAM~nlI;_1z1N5(B3}=fc3NhF>8|jgYFez zT`fQj?xnOTT!8hp0QnWAa1w6|u+A1B{jnE^{Fw@|-WDR!`3YQzb+-`F2ls$>A=ckQ zl;vms-RE{8*5N{QyZ-c+-q=E{$AyUA+P&mLtjmRHzPU4nSf7i~Ycuh;T_QzTr;AWn z_{ycbRg17*7a^0^?5Yp07h&BlLfd&K?|lp{!unl=rsaD}ch(hQ9WO%V_H8rnqeWQH zi;?HcN|}z;#aP#ik>=-&tJ@TcvA!3h)q;_6u2+k(&KIL@HQ%l`fyG$wi&2uB(;CUL zVyyech-S_-%PPkDUyKfXr#of|_JI<_x>>DW{6GozgA&xNe~{(2aS8T?5=5sJigyY2 zhZ1zQc8&d}yb|maC1~?6?(`d7CD<=Y&{Xi8Z69kX_Ki~1Av59YBUp<4qZH9cbDwl6 z_K{K)b@_8+&!JN6C#C2lHnb zlw!XrMS_<3XXFw~vG0_kKbKdE{ihUl#l#0bY%Rq;REnyJ@J!i(QtU^iD6u1m+&Nc@ zeW?u5*nHS?Ng4L1GGx8PgL~KJGVD`ji0-Bz_myG4Dnt8}>=f>+mto&3L+iV*p5J+~ z4Et9ZqIK?__YP&)$I4LuN`q(akIJy0m7&HAZ8|bzUn@gt$5mgo=agZ8D?|PRhiM6* z4EtOeqPz39Pi5He%Fr1h8Oe7OW!U%1k$nCZWmm3p?0@BGleWkuDWP)ggXM^8M9}rk z&T{OB<%rhUaA&1*?2F}yvXQHJ>zQ)wkL76IUMf&I7w(cr~G^hpKwh)V3^mB_-DUeT4<&nwZqUC32pU#~=RiC`T4fw`a(opPSjw-~L&d{Bw#=~s}w3UfjgnvcJ$t1vHA zAvP_#U{zsms6t=9w$q5G3iCr1njc-gD$Eg8C_|)y>&n$C%o9~;et2D~FjrKeTYCA2 z%L1z~UsR!U<7tLwNmZCLs?cHY_}k0MsxWU@f_|N%sVy6{EeWQ;rbfPJvB&uUa)E~|I{EE zKCbAjqcxa=YS5~9#`3cUHJFEL(A=Tv544iTTvUTTeV`|*8q7yEsLFw^ST&fFYS15i zufe=jgXU?vss?jY4YE>hxu*ZF2J=%5I@weGknLLy=BOI9?`|;-5o$3{)gqzQ+1003 z)?%)zMazWs}?nUdA3&XWG&{cT9kUvZ!fc1E#|ITi5U`PZ4SI?Q`@Xl#Cf>oE7#p*Jmb zhp)r@SBL2IIQ^jxb6_0`W1Tn`IaY^xunxI~4LPcE)MGBJN2c=&h4q*Z>k-|J!nfCB zPOL|?#e~*lN9r*z)}#5je5xLEV?COVDd!8jR*(6y9(^F_hE$I^vL2Opq$!OD)nlHl zM=^Hs7hfdTW3H@6f9$9p^JP7{mK+p1`l=ptW<8o;4f;}#d9xnv?xquK1Ln>KwE8yP zks2_6HXw#ocC<*=fH|}Q&EwU<2F#-ksP?$A>Nm{>%%u$|Y2eJ|;7bjdPaDu9jY1W+Hcc{Ik*vRrkArvBj({oG{5@zyb*J8BN_>Oy=Px~ zBj)2q)ViaEHia88CpRKmShen1ecT=OQ(*-hwOl*{SPyG@w4n@};j zY47RZgt@y3(bKO?LKEiiCgg5=IsHR%6Xx(HG{2liZNfa>gwzuChIaHfVJ>e%GDkEx z-_JH-K5s@VKPhc;U)+p2y&2J9tIHjCDMZ!6Wgp4X)>2foZab6;$SC2I>o#Q3r+(br72FWSqmuhz4V{6*w}^V`L<~WBsx- z_sKYyk;q)WlMk%AuZC@JWi zr{MfZLBHNs&~717aE_#)`Q!J)N*AIB!x=!u$^f=S~WuuU7odq~QEXLEZ`xH780aIEPXY?PkWS_&N&CqZDM_ z8T|Pdm4b6A1@$D;7oaIPpHk2~F7{J!PNkqaC0^PsGYZbD6tr*oeVQawaBig{I=#lR zQ*nN!BB6P~qT(D&MZULZk7}-=;yg=5EGqOn2vnSFsp$HuapT!-RGe?AsK56bEiqDY z&ZVMLnRLaX;=D^m^qYe})u=f4QW5>`mHH_v&c9TYy5?+3v=J5OU@8(Gd{;Dcm5TE) z75V3!q}?K<;#^Ec>~0jA#ZqxTrXu?4;j|YO=VU7SB2+(lJcx?(G8LWqs=yW%LB+Y5 zidyo^1g4XyI6qSnF&|i{I7d^_pXClJ&eK%1X&!&5I9F3q&`2t61%!(8H5JiUBO*Gf zIA>GQtp{;tlY>;8x2cF$%r{gfs5p01(YZ}wv|FN7oWENUWo$JKnOboUZ$-3l8GTHk z73c9*ls&(I-imX1D-xeKyjJu-_gqFPeuUwF_j4KlbN#`f{^uEt|Ly(p|9<}TAD?IX zub$samHN}#zt2(qt37|R(4K$Z1ONE^-rwg>mHy#*xXagU%`cucxoBu(di4_Rlz*NM zV_<-Pe)|7(@BiFi`2oPGHx%29Hd%WaBAO757#o`hmy8YE>W7&`xxjxId00M!qJD^Duy-gteEL z5IaO&t|yz{;s?Kp@u)&y0pQ&o#GB762ntR1D?3-L2RqIK*V1o^0Aq`9Q|qqHz$voW zGx*RJxOYX+KqE#P_NgXbJTAQhR>rjW`)A8QmiVOnCez*E>AOB_v*KR3|E$tQT;l+k z5gzSzhXB0v#WyNW5U}W}daF>;As90*yp$HA0E^CwY&hzp2+T=)8)lv#1&;o$Yu=R| zgI}8~JGRZKLb5?^>Jc$@2>$J%!fK!aEtk2CK7?w*#Eoe0&aYY!EjQqCO-UPWnd>K{ zr<{ZjgVB{4BD%of%=YY7lpdI{Q#)jKoCcfxbC=IH>BHF{F9R1}Hvqk{(FYRC&%u)H zlY^_K&cQ*EJt7mGhM+k7N{gF(9$XX)LoU`{fa83%i#cnIpnPaa^-SYM5I@AJrv1(s z+9yfnx4&J21MS`4+<8nPYQ>kwcH1w5rMT8I4FfZ%**{dW_K7*PR4wM#AYFmQzl_Oi zmtTW3G7q{A@2~)wpf?Il?=65K?t5>7mn8&#PtEy3Sb^RB{QEOAR&dQ{#nR$jYY13> z&;0OB8#twzp72=S7DCu+2BZ0|L!(7weD=rdpfKa*&6{lp&OM97Ui#QWoS<}fh?N64 ztr@m>t9t`%&bv9A9=-`imzX&{cin>h!>i=|w>pCR;U7!A#c#v!xNiBITkpUHe(swo zf_LFg_Dr=<1Z}Ie+bUwK6L_y(ZtRfn1XI&lJqgF0VNW;f;ysbP8`1P_vuM=*d91v<+HRA@(sVW=EQTM>LeU)B@raQ1b)4p<*#RC@o z$S-8c_W<4!6SK$G_d%Vf@0Im-Pv}%oew07v34V*~G%w{pfKhUB_bs=FP)e$2_;tby z9II3=eHD8I9HFT~w-~)4^{dIA6(7C9mhiXkI`dW(xFU%Z#_wH+vFBn~^?OBAL!1ywur<*FCz|ASmmTk&@P&(11zrEEDW_dK1 zOBnh?>a$tWn(zLwF7a8kplbjeiZ#=?B^U@9N|{U4GXo)H^Sf0+8bQ$Rwfo`8(I8lq zkm}_X5DeFD><4u4DWHVr-B$;m0+CdI%Q`3oR82hJxE^~3e%!L;D5mGo#2Yl&R`eVg zc~3_RI=+A@;b%^^5DLPs|NRY0dZ6qXaeic(c90~0(z4_JUC^$8cDV#SD z1w#*I^lik4<|e2HDfF=MZDnP@s5R=ioA;RLUF+9{B*Nm zMjS}>x~R)2$3yv|aJ`Jr@vu6LK{wbr0W8lyHSZTmgooTtsS&k_&>(x{ZkpM{?2d9sVJckH zlcUzIPlL@73JH?#X&`u_tabKV8Vrw%wC_Hh4qv%>6*_CvLC1C?Xi_c%q%1d`xs#Xy z(p`Pa9oA=peT@IH-oQ*CXMBD!us91U*s{b)kFvmWK(PBcUpBbhxgCANCmYl=J(Pbe z&w-s4wF$1FIbam_`5EJuT&VsX$!VFF3x2(#3cbpC@O9NdYxT!GK%WxH!uI*l;Sgs! zyR-m267TkFXB5B@Zl8PGPZmN@oxh{cY$5PJc@$(DQv}vQpUiu-iospF#z3951d1hJ z2>0ig0J1)ksk65fc0B61dLpa$L{^r@Zok_ra)0OoZ4ws@C|BU!!!kEhVhv z_MsMXpZMjQUaf=msb^&X0<}NUvg9RmNpO|dtVZdY6Hz1H(zoHy#lXw>kb${ zdj(t)CA*Z?yoSCP{Yk?wUPHnMfm7E+-+(cDn^0rY8~8rt^^r~PEvRic@35izE$k2d z*|1jU9pnl{WzG%01AQAM=Srt`U|@3hOkUjq%cET{z0K`_9&O%3%6jj?_yci>Y3@D9 z4^P}?j`{#c1!lIWAO8rW63tfze}07Vs}@d|Pj-U;W%t_1)J}-^sc}i?`2=EZHiLne zKY^p$ik;f^pCCR`?a7xdU2syRllE_`c0qt{j(Oa87dY>h7$=|ZhO^SaXT~eLp^E9n z={4IvL$_>}!m;O{q1^MFl`~HdY}k2ccZFjQI5mq7@_p|C`w)!}#>TxM9_Y%^^1c_k z<9*A{>h!^vca73(seO<%6IjWi(GOg{Ri{PC{jl4P$=OP40Az=L9UOZ-0I~9+3~>g7 zV9#5>CAe=8^846D$yQ&$=h<81;|xQPH>YT`JszV)aq zR~?3^MhPpH-eKU0x-G2a^bOkigKACrM?g%P%5dWB2wYVfaJY~<0`(2nV?j*c;o{AJ ztrObcfg?@fvwgyM$X&k7cO%0nsEOZ@8a_D+ESFo1T3(JqGxwsQI`%PGeWKyU2g5P= z6@8|Lw`dI9gTA@zt{MlKu349o>*Fw{sX8m!HV$)z3O-L|Cg7$?KrgHR1Qh7ygj)QV zfXO9a&$jFSfE{zqR!0hdK-;28rL^^vG##>xo^zvZ^~c874Nbx+Yw3bzT2r7iTyV21 ze+tM6{P}Jo(=cY>G5!1DG#EULSgAHS4eMUku_&IOf#fMYl@Cob;L`Tw=(4>(VKlno zmvd#G0y3u`5=rMC15TXMQ-bHN0)6O4$s}Hjg zf(|{AA$Dwp^W&o|`3Y>qQD*5g*Wa@dQUk)wE&S|6aosPCB7&WWj|@8q*6f62EW5Ez z6gzRVdf14mjrKa%GY+wG5H{tT<+*lo5ZlLXjSY=Bh~~M&afhF95YMhYdzw+kL8u6} zwH)^=Qkw#vJnt6^NkiIrf(#)ns& z5?{$pq!p^2nL5Z#M5HsY$X(cvSnNzum)WCW%&sJ^%Qpl@rOmP$Y zd-gt~tmYwpdOx%GKFCA-$f=D;ImbghY+f9_!o({;`xdG!)Ytz68H%g{o6Z& zXxp^O{&EA_Ry0?KQH{1=piMT>Dwdx(K*T=HisC1>ogZyF9L`S+Xmh+f^PHdX5k1P! z5X?^ma#gF2JfXFBOON=LNBqRW>+H?$?)*g2x6O4kcle2O&otdX+VK-N#Wk*IU*jj- zx{o{uV}3&0<-zlN`n2mRW3Xh?q`+ zh$C0tFUwfNM>HJpF&SauBdlXv;_eOc5^GA*+3l&k1b@c2^`&{d#7?bcT1P{93E|np za!Xx#iFJ!C7$r=2iQLBbM}v;>5|@|j7haI!B`l9FUGiuVtqz&q=qwxMA$~It+s{yV z2;=G7JqI#)h!Dvtn~%Oc#0{P&j2|s|h@Q+Bm-nml5K2;WJD8<-h_?62TlaGF5F_rd zTRQu=3Dwx6k0vX)iOMZPS8j)L6K zaVaWZIIWY5D2V5+{+P=}EbWsy@Aa6As9JNgu;3CGp{aI>&Fml-A*2;<86v<%eDocY zRvqLdDjpuQzFWdch^gma+2zMcq&Qk@n454CQt$26H}2;ok`G^0KhMue>~|6?km%wd zE{2LIxnyz>N8Z{?o8RLg26{6?dvrL6%ZqJlzKL-VL7W%bJg3^6iRRwGZPmB zWu&&aGZRsF_&8S3{6ofUh}1~)6W#M0eN{d&5qtAJ4Fcksh_&`#loBkN2nRmlWVby` zgrrjI{$q?xM9r0+OH1k)iP5lb$@a&L#EujJuV`IHB53I=7t8gGgp29*MW~;FVB%BQ z&X>wSNUrn!{K1ZaNV9z@mWCJzspKy|UD+83JnJ<0=tF@ccd8ywZ>3WD)vG^Cr zK6r1JeEb)P9%io$_&5hfE61vsEa$-a6iMSd(;TpFYuI<*e-4{6xsWE2IO>=HDw%T;8{xD>lD5jXy?)Qdl@wi5gm2Uo%T&bPRrU= zh&lx~yb}X@&89&4rB;I->lCd1!o6%XbP}djIvut5O+v@clf~h$e}LlD&Ihk-en4!d zT=3YEAF$qV?enR$3D~3Kd#p!e0&=c{Uf$?9NQ*Syu<;v*!&f4x#O`qz?38BQ(>?~` z#>1rJcgMhRbbFukrZLbssnFlrGz!NSU$#GdeH5HU-|Qm`jlw2M$*UXn-$DPX+E=*# z9kvPti-?GPhcAW6Yd4TbKt{i7$NRyf?4;K3d?;VTp zFg_lJ6C06p!t_^Qv>4;{yZ99bKkI!hqz-|o(d$mp!$ZKPQMRuk`wN^Bh?7s*@&)2l z_Ut7^48ld9Nn^=1gK+qa7H@Un0C@G}SOzZ}0IERHn;8FoP^v`tzOLwpvh78zPS5(_ zwoBuuO(K1e@ThcgOJXmG9Tom+FWU>2Jy#k;Dth4Q3L%yjjUJfT@NC$j>odrHkmwSz z{S0@uFP3fO{0y8n?u1oXH{6wfb$Z#JZV>%`@>3?Y3zXMnNvKKHu9WbplDl%!&4&qxpcsSp` zgO8PcGN)AD!4k2E$d2N-AQZ7}AD_%ysIn`*+no3YHX1AmS+nU4gl+qnIvMdAHs7e| zHxPae=YA=fYDB$)Q2Ecg&0?=$YWev4)|5676)Ib2w5JVroOUE`R<*+I)I;X0PPBri z)rzVEJyalh_E`4WQ{fU39JGv|3dM%(MY*vQcxz&I_LMvYJeLTKef~g(C9eb9bR5YL z_Ez2~eJvSGRyoa%m5?BXC+_EJBNEJsQi{d7NDyqX^~fHh7RX4N=IYOBhPmb$Hv|4= z$cgxA5puZ+x@4o=p4TA89+3ec6hez_i8nm_lu zzN>@o!0=0_RqMcHaZ}M$c`bAs>Hj*tuNH2u+qB;`y9PXje&4r{sevQA17?J>s^R-< z$1oqcYFHU4VQqZ~IXr*|k1R}bse}z$tVUqTrVGM04hFy^1E$&eS3>-;` z+Wm#Je$(o@rcwwqt)oaivj7ZkFL|&^r~nS0Jh`plVLmV&z00#}E)OR2SJt>$v|U9~zvKo&;x~ z2gf@LCc#Z!$7@nKG`Lt~qjAH3qhP5wx**9SuTRVs)pC zqJiABGiG`$3W{QjOmf|#z=x{I?Y1fkw3H3x&!t4dp8k)9M-?L>YdF6MeTsmN@2&Cj zju9~F+~lCKIs(3jDmj@Igu|3c_tM8_!lBVf#;uq!9JV`ViT#KT1I9Zo8nWtPu+A(a zc6K}zK4xqTeX1G?6`yVotx0$Rt{e|86>z?QgQ2%>9yfmuHFs?6 z(Q1ul&CwwszE>$mP9Ox*xu)Cv-JXK|&kntw>0s~}e7kp(WiU`X8=Z8%1i`&ES9aT* z1c5S_sgZksAe8Rfwb<1x5Wai&mYaPKfO_Z7eBK)Y@SWw2Zx>4dB-{)z{^;WmM|@2* zSvL5?laYx?^>jaA^ly!;Q1An$s~#`)+n<2!m^0IA>nD&{EIc~L^8`L^%&2FH_l0PR z*;^?^UuZaOVRp062S(nvG!5MMfwqyx0}fI?u+)cfLa6OAW3uG=f^6@xd7O3x}A8$5#IwgDEQ5-*SzTlw|3uopOFYX4&QdCmipp$<->n4`s(9+&9bKhoF;=VWsIFFnS@t zENi<5Y{}w3&6eQ~)!r<((q(DehjY8~3h#lZ?YhoY3im+6Y(pNk#SKT5BaF%=Q_cx@H&r0r=38Y*;RY%@4FEF(L$6r`Ys43c^ajw+=VKgey-g^ zcfh*$$yaarJMgl4t@_nxw}Fs1Ure5M1lBtT8FEfJ!fBbv*UL+8!A7&@qQp(N;LJ

k7Cz|TRk8eH z148K!9j6^^;PM@{4Z&PC(5rpOTkW|u+%S8;Zp{vBSaH;D^RHSfICC+OVZzV~UWZ39 zGR#?m4u#!$-BU|YEo*wEw#O2JNljUg-&w#Lw(PhXM+>+w4$ng^mKZCTLt&ZKXtkqrl3ZpG zYcE=1X?Yn083Z(j+fBjc%`yS&qo(lHaz-L1&jf}#U-Q{)HGyD#ldSjQmmpYJdwbWa zOAzqx%cHw~#-MpOgE@1tF_;GHdEW842%R=!yStWNgy{RFYukg2fZ5+<+<3haxZAh# zm?c~Quj`HLVq`BsPn3Sg&&u=gaLw-xSGCWBRmphY&ptyqTQjvW>Xsq&-Z7JL zYnrV`py$b~g{bfm+Is2ITidM^;OUu%Lmt(K!R5G~&E0i}p{kQ%-SQiU;LR~viM?;- zVdj(c_MiLYVcYTeAL}CsnBFC)5=L7-32jmu(62WHouWeZbxgNjp3 zu5gVk{c1rcH>z{s_st6kWwo@?YRTKIoQ8GP_`Xv zmxdlyp4$duGmEy}Ik*kjd)=E-Ii=yDk(S@4FiD6b$39;LTVbaJLMa0hAfS1UgZxY! zTpw8^8taL{bo*gLGoj6(@*>597~Tk)p;xcc=D7*2o5MN(m1+yE$IIEMSXf?HB@rQhwZw%3Vz>zG4ZKx1x$?& zQ+h;}L&C{q@dJ+pfSf23mBz6Y*6Fyg%=j(_nMDs#^43Mr%u;zXPv!tC>1NAy`fPAO!c_fZ9CwH`fwpg4<2~6(W1$IBc^dNw8n(FjB)oqR5sD`?4<6`6Z5%BZtXaB{>8RA; z#LW&=v(NCYG_?)!{%qsq;3J`ovnK_ebgB_0&QqW#r4W56jGeWh{lD~K2Zc^sN<2ps V;%U!+ZoZ95rkjuN5dC!U-vFWH8jt`0 literal 0 HcmV?d00001 diff --git a/test/test_flow1d_coupling.py b/test/test_flow1d_coupling.py new file mode 100755 index 000000000..d0a593f55 --- /dev/null +++ b/test/test_flow1d_coupling.py @@ -0,0 +1,66 @@ +#!/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' + +# Save output with hdf5 +with h5py.File("flow1d_coupling.sol.h5", "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 From cf3a272c3a836757e2a3622602cc19247d433e59 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Fri, 23 Aug 2024 13:41:58 -0500 Subject: [PATCH 7/8] Specify version of mpi4py Latest doesn't seem to want to build in this environment. --- docker/test/Dockerfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docker/test/Dockerfile b/docker/test/Dockerfile index 3ec07450a..50e665964 100644 --- a/docker/test/Dockerfile +++ b/docker/test/Dockerfile @@ -151,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 From a1306c2360d34752c5bd473336101038953b455d Mon Sep 17 00:00:00 2001 From: Joshua Bigman Date: Fri, 23 Aug 2024 17:52:46 -0500 Subject: [PATCH 8/8] Minor changes to fix regression test failure. --- .gitattributes | 2 +- test/Makefile.am | 2 +- test/test_flow1d_coupling.py | 3 ++- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/.gitattributes b/.gitattributes index 73e0aca5d..e99208be6 100644 --- a/.gitattributes +++ b/.gitattributes @@ -27,4 +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.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/test/Makefile.am b/test/Makefile.am index 5d0ef1593..9ee1c876e 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -35,7 +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.sol.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 \ diff --git a/test/test_flow1d_coupling.py b/test/test_flow1d_coupling.py index d0a593f55..36ac916b8 100755 --- a/test/test_flow1d_coupling.py +++ b/test/test_flow1d_coupling.py @@ -58,8 +58,9 @@ 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("flow1d_coupling.sol.h5", "w") as f: +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)