From 6af2ee8099dbed5afd0a9681a818b1292d5aa0b8 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Wed, 29 Jan 2025 16:10:02 -0600 Subject: [PATCH 01/32] Preliminary implementation of bte coupling for low-mach --- src/M2ulPhyS2Boltzmann.cpp | 1 - src/reactingFlow.cpp | 60 ++++++++++++++++++++++++++++++++++++++ src/reactingFlow.hpp | 6 ++++ 3 files changed, 66 insertions(+), 1 deletion(-) diff --git a/src/M2ulPhyS2Boltzmann.cpp b/src/M2ulPhyS2Boltzmann.cpp index 8bdf6a5cc..2286035e9 100644 --- a/src/M2ulPhyS2Boltzmann.cpp +++ b/src/M2ulPhyS2Boltzmann.cpp @@ -36,7 +36,6 @@ #include "M2ulPhyS.hpp" #include "tps2Boltzmann.hpp" -// CPU version (just for starting up) void M2ulPhyS::push(TPS::Tps2Boltzmann &interface) { assert(interface.IsInitialized()); diff --git a/src/reactingFlow.cpp b/src/reactingFlow.cpp index 9e37af02c..2a4c49480 100644 --- a/src/reactingFlow.cpp +++ b/src/reactingFlow.cpp @@ -41,6 +41,7 @@ #include "loMach.hpp" #include "loMach_options.hpp" #include "radiation.hpp" +#include "tps2Boltzmann.hpp" using namespace mfem; using namespace mfem::common; @@ -363,6 +364,11 @@ ReactingFlow::ReactingFlow(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, tem tpsP_->getRequiredInput((basepath + "/radius").c_str(), R); rxnModelParamsHost.push_back(Vector({R})); + } else if (model == "bte") { + config.reactionModels[r - 1] = GRIDFUNCTION_RXN; + int index; + tpsP->getRequiredInput((basepath + "/bte/index").c_str(), index); + config.chemistryInput.reactionInputs[r - 1].indexInput = index; } else { grvy_printf(GRVY_ERROR, "\nUnknown reaction_model -> %s", model.c_str()); exit(ERROR); @@ -2835,6 +2841,60 @@ double species_uniform(const Vector &coords, double t) { return yn; } +void ReactingFlow::push(TPS::Tps2Boltzmann &interface) { + assert(interface.IsInitialized()); + + const int nscalardofs(vfes->GetNDofs()); + + mfem::ParGridFunction *species = + new mfem::ParGridFunction(&interface.NativeFes(TPS::Tps2Boltzmann::Index::SpeciesDensities)); + + double *species_data = species->HostWrite(); + + const double *dataRho = rn_.HostRead(); + const double *dataY = Yn_.HostRead(); + + double state_local[gpudata::MAXEQUATIONS]; + double species_local[gpudata::MAXSPECIES]; + + for (int i = 0; i < nscalardofs; i++) { + for (int asp = 0; asp < nActiveSpecies_; asp++) + state_local[asp] =dataRho[i]*dataY[i+asp*nscalardofs]; + + mixture_->computeNumberDensities(state_local, species_local); + + + for (int sp = 0; sp < interface.Nspecies(); sp++) + species_data[i + sp * nscalardofs] = AVOGADRONUMBER * species_local[sp]; + } + + interface.interpolateFromNativeFES(*species, TPS::Tps2Boltzmann::Index::SpeciesDensities); + interface.interpolateFromNativeFES(Tn_gf_, TPS::Tps2Boltzmann::Index::HeavyTemperature); + interface.interpolateFromNativeFES(Tn_gf_, TPS::Tps2Boltzmann::Index::ElectronTemperature); + + interface.setTimeStep(this->dt_); + interface.setCurrentTime(this->time_); + + delete species; + delete heavyTemperature; + delete electronTemperature; +} + +void ReactingFlow::fetch(TPS::Tps2Boltzmann &interface) { + mfem::ParFiniteElementSpace *reaction_rates_fes(&(interface.NativeFes(TPS::Tps2Boltzmann::Index::ReactionRates))); + externalReactionRates_.reset(new mfem::ParGridFunction(reaction_rates_fes)); + interface.interpolateToNativeFES(*externalReactionRates_, TPS::Tps2Boltzmann::Index::ReactionRates); +#if defined(_CUDA_) || defined(_HIP_) + const double *data(externalReactionRates_->Read()); + int size(externalReactionRates_->FESpace()->GetNDofs()); + assert(externalReactionRates_->FESpace()->GetOrdering() == mfem::Ordering::byNODES); + gpu::deviceSetChemistryReactionData<<<1, 1>>>(data, size, chemistry_); +#else + chemistry_->setGridFunctionRates(*externalReactionRates_); +#endif +} + + #if 0 void ReactingFlow::uniformInlet() { int myRank; diff --git a/src/reactingFlow.hpp b/src/reactingFlow.hpp index 1387ffced..4e19ac325 100644 --- a/src/reactingFlow.hpp +++ b/src/reactingFlow.hpp @@ -36,6 +36,7 @@ // forward-declaration for Tps support class namespace TPS { class Tps; +class Tps2Boltzmann; } #include @@ -102,6 +103,8 @@ class ReactingFlow : public ThermoChemModelBase { PerfectMixture *mixture_ = NULL; ArgonMixtureTransport *transport_ = NULL; Chemistry *chemistry_ = NULL; + // External reaction rates when chemistry is implemented using the BTE option + std::unique_ptr externalReactionRates_; std::vector speciesNames_; std::map atomMap_; @@ -415,5 +418,8 @@ class ReactingFlow : public ThermoChemModelBase { void evalSubstepNumber(); void readTableWrapper(std::string inputPath, TableInput &result); + + void push(TPS::Tps2Boltzmann &interface); + void fetch(TPS::Tps2Boltzmann &interface); }; #endif // REACTINGFLOW_HPP_ From 6b37d97891dcb5cdbddd2388fa40386279ce06a4 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Wed, 29 Jan 2025 16:26:30 -0600 Subject: [PATCH 02/32] Fix compilation errors --- src/reactingFlow.cpp | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/reactingFlow.cpp b/src/reactingFlow.cpp index 2a4c49480..1969db986 100644 --- a/src/reactingFlow.cpp +++ b/src/reactingFlow.cpp @@ -365,10 +365,10 @@ ReactingFlow::ReactingFlow(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, tem rxnModelParamsHost.push_back(Vector({R})); } else if (model == "bte") { - config.reactionModels[r - 1] = GRIDFUNCTION_RXN; + reactionModels[r - 1] = GRIDFUNCTION_RXN; int index; - tpsP->getRequiredInput((basepath + "/bte/index").c_str(), index); - config.chemistryInput.reactionInputs[r - 1].indexInput = index; + tpsP_->getRequiredInput((basepath + "/bte/index").c_str(), index); + chemistryInput_.reactionInputs[r - 1].indexInput = index; } else { grvy_printf(GRVY_ERROR, "\nUnknown reaction_model -> %s", model.c_str()); exit(ERROR); @@ -2844,7 +2844,7 @@ double species_uniform(const Vector &coords, double t) { void ReactingFlow::push(TPS::Tps2Boltzmann &interface) { assert(interface.IsInitialized()); - const int nscalardofs(vfes->GetNDofs()); + const int nscalardofs(vfes_->GetNDofs()); mfem::ParGridFunction *species = new mfem::ParGridFunction(&interface.NativeFes(TPS::Tps2Boltzmann::Index::SpeciesDensities)); @@ -2859,7 +2859,7 @@ void ReactingFlow::push(TPS::Tps2Boltzmann &interface) { for (int i = 0; i < nscalardofs; i++) { for (int asp = 0; asp < nActiveSpecies_; asp++) - state_local[asp] =dataRho[i]*dataY[i+asp*nscalardofs]; + state_local[asp] = dataRho[i]*dataY[i+asp*nscalardofs]; mixture_->computeNumberDensities(state_local, species_local); @@ -2876,8 +2876,6 @@ void ReactingFlow::push(TPS::Tps2Boltzmann &interface) { interface.setCurrentTime(this->time_); delete species; - delete heavyTemperature; - delete electronTemperature; } void ReactingFlow::fetch(TPS::Tps2Boltzmann &interface) { From b9fb3507156206ff06d70b6af64d1c7d0a1f8b34 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Tue, 18 Feb 2025 20:25:39 -0600 Subject: [PATCH 03/32] Working one-way --- src/cycle_avg_joule_coupling.cpp | 36 ++++++++++++++++++++++++++++++-- src/loMach.hpp | 4 ++++ src/reactingFlow.cpp | 2 +- src/reactingFlow.hpp | 4 ++-- src/thermo_chem_base.hpp | 11 ++++++++++ src/tps2Boltzmann.cpp | 9 ++++++-- src/tps2Boltzmann.hpp | 1 + 7 files changed, 60 insertions(+), 7 deletions(-) diff --git a/src/cycle_avg_joule_coupling.cpp b/src/cycle_avg_joule_coupling.cpp index 07ab637a1..4456a72c6 100644 --- a/src/cycle_avg_joule_coupling.cpp +++ b/src/cycle_avg_joule_coupling.cpp @@ -348,7 +348,23 @@ void CycleAvgJouleCoupling::interpElectricFieldFromEMToFlow() { if (flow_fespace->IsDGSpace()) { efieldR_->SetFromTrueDofs(interp_vals); } else { - assert(false); + Array vdofs; + Vector elem_dof_vals; + int n0 = 0; + const int NE = flow_solver_->getMesh()->GetNE(); + for (int i = 0; i < NE; i++) { + flow_fespace->GetElementDofs(i, vdofs); + const int nsp = flow_fespace->GetFE(i)->GetNodes().GetNPoints(); + assert(nsp == vdofs.Size()); + elem_dof_vals.SetSize(nsp); + for (int j = 0; j < nsp; j++) { + elem_dof_vals(j) = interp_vals(n0 + j); + } + efieldR_->SetSubVector(vdofs, elem_dof_vals); + n0 += nsp; + } + efieldR_->SetTrueVector(); + efieldR_->SetFromTrueVector(); } efieldR_->HostRead(); @@ -357,7 +373,23 @@ void CycleAvgJouleCoupling::interpElectricFieldFromEMToFlow() { if (flow_fespace->IsDGSpace()) { efieldI_->SetFromTrueDofs(interp_vals); } else { - assert(false); + Array vdofs; + Vector elem_dof_vals; + int n0 = 0; + const int NE = flow_solver_->getMesh()->GetNE(); + for (int i = 0; i < NE; i++) { + flow_fespace->GetElementDofs(i, vdofs); + const int nsp = flow_fespace->GetFE(i)->GetNodes().GetNPoints(); + assert(nsp == vdofs.Size()); + elem_dof_vals.SetSize(nsp); + for (int j = 0; j < nsp; j++) { + elem_dof_vals(j) = interp_vals(n0 + j); + } + efieldI_->SetSubVector(vdofs, elem_dof_vals); + n0 += nsp; + } + efieldI_->SetTrueVector(); + efieldI_->SetFromTrueVector(); } efieldI_->HostRead(); #else diff --git a/src/loMach.hpp b/src/loMach.hpp index ba5a9066a..dcff8254e 100644 --- a/src/loMach.hpp +++ b/src/loMach.hpp @@ -57,6 +57,7 @@ // forward-declaration for Tps support class namespace TPS { class Tps; +class Tps2Boltzmann; } #include @@ -258,6 +259,9 @@ class LoMachSolver : public TPS::PlasmaSolver { void evaluatePlasmaConductivityGF() override { thermo_->evaluatePlasmaConductivityGF(); } mfem::ParGridFunction *getJouleHeatingGF() override { return thermo_->getJouleHeatingGF(); } + + void push(TPS::Tps2Boltzmann &interface) override { thermo_->push(interface); } + void fetch(TPS::Tps2Boltzmann &interface) override { thermo_->fetch(interface); } }; #endif // LOMACH_HPP_ diff --git a/src/reactingFlow.cpp b/src/reactingFlow.cpp index 1969db986..93d9952d1 100644 --- a/src/reactingFlow.cpp +++ b/src/reactingFlow.cpp @@ -427,7 +427,7 @@ ReactingFlow::ReactingFlow(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, tem equilibriumConstantParams[p + r * gpudata::MAXCHEMPARAMS]; } - if (reactionModels[r] != TABULATED_RXN) { + if (reactionModels[r] != TABULATED_RXN && reactionModels[r] != GRIDFUNCTION_RXN) { assert(rxn_param_idx < rxnModelParamsHost.size()); chemistryInput_.reactionInputs[r].modelParams = rxnModelParamsHost[rxn_param_idx].Read(); rxn_param_idx += 1; diff --git a/src/reactingFlow.hpp b/src/reactingFlow.hpp index 4e19ac325..a4219bb50 100644 --- a/src/reactingFlow.hpp +++ b/src/reactingFlow.hpp @@ -419,7 +419,7 @@ class ReactingFlow : public ThermoChemModelBase { void evalSubstepNumber(); void readTableWrapper(std::string inputPath, TableInput &result); - void push(TPS::Tps2Boltzmann &interface); - void fetch(TPS::Tps2Boltzmann &interface); + void push(TPS::Tps2Boltzmann &interface) override; + void fetch(TPS::Tps2Boltzmann &interface) override; }; #endif // REACTINGFLOW_HPP_ diff --git a/src/thermo_chem_base.hpp b/src/thermo_chem_base.hpp index 78477259c..003934964 100644 --- a/src/thermo_chem_base.hpp +++ b/src/thermo_chem_base.hpp @@ -39,6 +39,7 @@ namespace TPS { class Tps; +class Tps2Boltzmann; } class IODataOrganizer; @@ -197,6 +198,16 @@ class ThermoChemModelBase { std::cout << "ERROR: " << __func__ << " remains unimplemented" << std::endl; exit(1); } + + virtual void push(TPS::Tps2Boltzmann &interface) { + std::cout << "ERROR: " << __func__ << " remains unimplemented" << std::endl; + exit(1); + } + + virtual void fetch(TPS::Tps2Boltzmann &interface) { + std::cout << "ERROR: " << __func__ << " remains unimplemented" << std::endl; + exit(1); + } }; /** diff --git a/src/tps2Boltzmann.cpp b/src/tps2Boltzmann.cpp index 689fded06..402721958 100644 --- a/src/tps2Boltzmann.cpp +++ b/src/tps2Boltzmann.cpp @@ -79,7 +79,7 @@ void idenity_fun(const Vector &x, Vector &out) { } Tps2Boltzmann::Tps2Boltzmann(Tps *tps) - : NIndexes(7), tps_(tps), all_fes_(nullptr), save_to_paraview_dc(false), paraview_dc(nullptr) { + : NIndexes(7), tps_(tps), all_fes_(nullptr), use_h1fec_(false), save_to_paraview_dc(false), paraview_dc(nullptr) { // Assert we have a couple solver; assert(tps->isFlowEMCoupled()); @@ -92,6 +92,7 @@ Tps2Boltzmann::Tps2Boltzmann(Tps *tps) tps->getRequiredInput("em/current_frequency", EfieldAngularFreq_); EfieldAngularFreq_ *= 2. * M_PI; + use_h1fec_ = tps->getInput("boltzmannInterface/h1fec", false); save_to_paraview_dc = tps->getInput("boltzmannInterface/save_to_paraview", false); offsets.SetSize(NIndexes + 1); @@ -119,7 +120,11 @@ int Tps2Boltzmann::_countBTEReactions() { void Tps2Boltzmann::init(TPS::PlasmaSolver *flowSolver) { std::cout << "Tps2Boltzmann::init is called" << std::endl; mfem::ParMesh *pmesh(flowSolver->getMesh()); - fec_ = new mfem::L2_FECollection(order_, pmesh->Dimension(), basis_type_); + if( use_h1fec_) + fec_ = new mfem::H1_FECollection(order_, pmesh->Dimension(), basis_type_); + else + fec_ = new mfem::L2_FECollection(order_, pmesh->Dimension(), basis_type_); + switch (pmesh->Dimension()) { case 2: nEfieldComps_ = 2; diff --git a/src/tps2Boltzmann.hpp b/src/tps2Boltzmann.hpp index 5dc39adb6..2d0c73bc0 100644 --- a/src/tps2Boltzmann.hpp +++ b/src/tps2Boltzmann.hpp @@ -170,6 +170,7 @@ class Tps2Boltzmann { double timestep_; double currentTime_; + bool use_h1fec_; bool save_to_paraview_dc; mfem::ParaViewDataCollection *paraview_dc; std::vector reaction_eqs_; From fb70d2af605bd689139266c07ab0471d727b99e8 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Tue, 18 Feb 2025 20:55:48 -0600 Subject: [PATCH 04/32] Fix asyle --- src/tps2Boltzmann.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tps2Boltzmann.cpp b/src/tps2Boltzmann.cpp index 402721958..8de6dc79d 100644 --- a/src/tps2Boltzmann.cpp +++ b/src/tps2Boltzmann.cpp @@ -120,7 +120,7 @@ int Tps2Boltzmann::_countBTEReactions() { void Tps2Boltzmann::init(TPS::PlasmaSolver *flowSolver) { std::cout << "Tps2Boltzmann::init is called" << std::endl; mfem::ParMesh *pmesh(flowSolver->getMesh()); - if( use_h1fec_) + if (use_h1fec_) fec_ = new mfem::H1_FECollection(order_, pmesh->Dimension(), basis_type_); else fec_ = new mfem::L2_FECollection(order_, pmesh->Dimension(), basis_type_); From b2fdca6c9342dae946befcf4b90c636b08fb52bc Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Wed, 19 Feb 2025 10:11:37 -0600 Subject: [PATCH 05/32] Bug fix in density computations --- src/reactingFlow.cpp | 5 ++--- src/tps-time-loop.py | 34 +++++++++++++++++++++++++++++++++- 2 files changed, 35 insertions(+), 4 deletions(-) diff --git a/src/reactingFlow.cpp b/src/reactingFlow.cpp index 93d9952d1..565ae4d8c 100644 --- a/src/reactingFlow.cpp +++ b/src/reactingFlow.cpp @@ -2858,12 +2858,11 @@ void ReactingFlow::push(TPS::Tps2Boltzmann &interface) { double species_local[gpudata::MAXSPECIES]; for (int i = 0; i < nscalardofs; i++) { + state_local[0] = dataRho[i]; for (int asp = 0; asp < nActiveSpecies_; asp++) - state_local[asp] = dataRho[i]*dataY[i+asp*nscalardofs]; - + state_local[dim_ + 2 + asp] = dataRho[i]*dataY[i+asp*nscalardofs]; mixture_->computeNumberDensities(state_local, species_local); - for (int sp = 0; sp < interface.Nspecies(); sp++) species_data[i + sp * nscalardofs] = AVOGADRONUMBER * species_local[sp]; } diff --git a/src/tps-time-loop.py b/src/tps-time-loop.py index af2efc6c8..347ac924f 100755 --- a/src/tps-time-loop.py +++ b/src/tps-time-loop.py @@ -5,6 +5,38 @@ from mpi4py import MPI +class NullSolver: + def __init__(self): + self.species_densities = None + self.efield = None + self.heavy_temperature = None + #Reaction 1: 'Ar + E => Ar.+1 + 2 E', + #Reaction 2: 'Ar.+1 + 2 E => Ar + E' + + + def fetch(self, interface): + n_reactions =interface.nComponents(libtps.t2bIndex.ReactionRates) + for r in range(n_reactions): + print("Reaction ", r+1, ": ", interface.getReactionEquation(r)) + self.species_densities = np.array(interface.HostRead(libtps.t2bIndex.SpeciesDensities), copy=False) + self.efield = np.array(interface.HostRead(libtps.t2bIndex.ElectricField), copy=False) + self.heavy_temperature = np.array(interface.HostRead(libtps.t2bIndex.HeavyTemperature), copy=False) + + efieldAngularFreq = interface.EfieldAngularFreq() + print("Species densities:", self.species_densities.min(), " ", self.species_densities.max()) + print("Heavy Temp:", self.heavy_temperature.min(), " ", self.heavy_temperature.max()) + print("Efield:", self.efield.min(), " ", self.efield.max()) + print("Electric field angular frequency: ", efieldAngularFreq) + + + + def solve(self): + pass + + def push(self, interface): + rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False) + rates[:] = 0. + class ArrheniusSolver: def __init__(self): self.UNIVERSALGASCONSTANT = 8.3144598; # J * mol^(-1) * K^(-1) @@ -61,7 +93,7 @@ def push(self, interface): tps.chooseSolver() tps.initialize() -boltzmann = ArrheniusSolver() +boltzmann = NullSolver() interface = libtps.Tps2Boltzmann(tps) tps.initInterface(interface) From 3ff62d8915f7ed6b6f18c1b631666dcb0dcda267 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Wed, 19 Feb 2025 14:50:09 -0600 Subject: [PATCH 06/32] Fix bug in parallel --- src/reactingFlow.cpp | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/src/reactingFlow.cpp b/src/reactingFlow.cpp index 565ae4d8c..10027316a 100644 --- a/src/reactingFlow.cpp +++ b/src/reactingFlow.cpp @@ -2844,12 +2844,15 @@ double species_uniform(const Vector &coords, double t) { void ReactingFlow::push(TPS::Tps2Boltzmann &interface) { assert(interface.IsInitialized()); - const int nscalardofs(vfes_->GetNDofs()); + //const int nscalardofs(vfes_->GetNDofs()); mfem::ParGridFunction *species = new mfem::ParGridFunction(&interface.NativeFes(TPS::Tps2Boltzmann::Index::SpeciesDensities)); - double *species_data = species->HostWrite(); + std::cout << sDofInt_ << " " << interface.Nspecies() << std::endl << std::flush; + MPI_Barrier(pmesh_->GetComm()); + mfem::Vector speciesInt(sDofInt_*interface.Nspecies()); + double *species_data = speciesInt.HostWrite(); const double *dataRho = rn_.HostRead(); const double *dataY = Yn_.HostRead(); @@ -2857,16 +2860,28 @@ void ReactingFlow::push(TPS::Tps2Boltzmann &interface) { double state_local[gpudata::MAXEQUATIONS]; double species_local[gpudata::MAXSPECIES]; - for (int i = 0; i < nscalardofs; i++) { + for (int i = 0; i < gpudata::MAXEQUATIONS; ++i) + state_local[i] = 0.; + + for (int i = 0; i < gpudata::MAXSPECIES; ++i) + species_local[i] = 0.; + + for (int i = 0; i < sDofInt_; i++) { state_local[0] = dataRho[i]; for (int asp = 0; asp < nActiveSpecies_; asp++) - state_local[dim_ + 2 + asp] = dataRho[i]*dataY[i+asp*nscalardofs]; + state_local[dim_ + 2 + asp] = dataRho[i]*dataY[i+asp*sDofInt_]; mixture_->computeNumberDensities(state_local, species_local); for (int sp = 0; sp < interface.Nspecies(); sp++) - species_data[i + sp * nscalardofs] = AVOGADRONUMBER * species_local[sp]; + species_data[i + sp * sDofInt_] = AVOGADRONUMBER * species_local[sp]; } + std::cout << "species->SetFromTrueDofs(speciesInt);" << std::endl << std::flush; + MPI_Barrier(pmesh_->GetComm()); + species->SetFromTrueDofs(speciesInt); + std::cout << "Tn_gf_.SetFromTrueDofs(Tn_);" << std::endl << std::flush; + MPI_Barrier(pmesh_->GetComm()); + Tn_gf_.SetFromTrueDofs(Tn_); interface.interpolateFromNativeFES(*species, TPS::Tps2Boltzmann::Index::SpeciesDensities); interface.interpolateFromNativeFES(Tn_gf_, TPS::Tps2Boltzmann::Index::HeavyTemperature); interface.interpolateFromNativeFES(Tn_gf_, TPS::Tps2Boltzmann::Index::ElectronTemperature); From 5dfbd0afe9c88eee393d9f9c6297fccdfe79df08 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Wed, 19 Feb 2025 14:51:29 -0600 Subject: [PATCH 07/32] Improve output to screen --- src/reactingFlow.cpp | 6 ------ src/tps-time-loop.py | 30 ++++++++++++++++++------------ src/tps2Boltzmann.cpp | 1 - 3 files changed, 18 insertions(+), 19 deletions(-) diff --git a/src/reactingFlow.cpp b/src/reactingFlow.cpp index 10027316a..015fcddc5 100644 --- a/src/reactingFlow.cpp +++ b/src/reactingFlow.cpp @@ -2849,8 +2849,6 @@ void ReactingFlow::push(TPS::Tps2Boltzmann &interface) { mfem::ParGridFunction *species = new mfem::ParGridFunction(&interface.NativeFes(TPS::Tps2Boltzmann::Index::SpeciesDensities)); - std::cout << sDofInt_ << " " << interface.Nspecies() << std::endl << std::flush; - MPI_Barrier(pmesh_->GetComm()); mfem::Vector speciesInt(sDofInt_*interface.Nspecies()); double *species_data = speciesInt.HostWrite(); @@ -2876,11 +2874,7 @@ void ReactingFlow::push(TPS::Tps2Boltzmann &interface) { species_data[i + sp * sDofInt_] = AVOGADRONUMBER * species_local[sp]; } - std::cout << "species->SetFromTrueDofs(speciesInt);" << std::endl << std::flush; - MPI_Barrier(pmesh_->GetComm()); species->SetFromTrueDofs(speciesInt); - std::cout << "Tn_gf_.SetFromTrueDofs(Tn_);" << std::endl << std::flush; - MPI_Barrier(pmesh_->GetComm()); Tn_gf_.SetFromTrueDofs(Tn_); interface.interpolateFromNativeFES(*species, TPS::Tps2Boltzmann::Index::SpeciesDensities); interface.interpolateFromNativeFES(Tn_gf_, TPS::Tps2Boltzmann::Index::HeavyTemperature); diff --git a/src/tps-time-loop.py b/src/tps-time-loop.py index 347ac924f..cd5018298 100755 --- a/src/tps-time-loop.py +++ b/src/tps-time-loop.py @@ -5,8 +5,13 @@ from mpi4py import MPI +def master_print(comm: MPI.Comm, *args, **kwargs) -> None: + if comm.rank == 0: + print(*args, **kwargs) + class NullSolver: - def __init__(self): + def __init__(self, comm): + self.comm = comm self.species_densities = None self.efield = None self.heavy_temperature = None @@ -17,16 +22,16 @@ def __init__(self): def fetch(self, interface): n_reactions =interface.nComponents(libtps.t2bIndex.ReactionRates) for r in range(n_reactions): - print("Reaction ", r+1, ": ", interface.getReactionEquation(r)) + master_print(self.comm, "Reaction ", r+1, ": ", interface.getReactionEquation(r)) self.species_densities = np.array(interface.HostRead(libtps.t2bIndex.SpeciesDensities), copy=False) self.efield = np.array(interface.HostRead(libtps.t2bIndex.ElectricField), copy=False) self.heavy_temperature = np.array(interface.HostRead(libtps.t2bIndex.HeavyTemperature), copy=False) efieldAngularFreq = interface.EfieldAngularFreq() - print("Species densities:", self.species_densities.min(), " ", self.species_densities.max()) - print("Heavy Temp:", self.heavy_temperature.min(), " ", self.heavy_temperature.max()) - print("Efield:", self.efield.min(), " ", self.efield.max()) - print("Electric field angular frequency: ", efieldAngularFreq) + master_print(self.comm,"Species densities:", self.species_densities.min(), " ", self.species_densities.max()) + master_print(self.comm,"Heavy Temp:", self.heavy_temperature.min(), " ", self.heavy_temperature.max()) + master_print(self.comm,"Efield:", self.efield.min(), " ", self.efield.max()) + master_print(self.comm,"Electric field angular frequency: ", efieldAngularFreq) @@ -38,7 +43,8 @@ def push(self, interface): rates[:] = 0. class ArrheniusSolver: - def __init__(self): + def __init__(self,comm): + self.comm = comm self.UNIVERSALGASCONSTANT = 8.3144598; # J * mol^(-1) * K^(-1) self.species_densities = None self.efield = None @@ -53,13 +59,13 @@ def __init__(self): def fetch(self, interface): n_reactions =interface.nComponents(libtps.t2bIndex.ReactionRates) for r in range(n_reactions): - print("Reaction ", r+1, ": ", interface.getReactionEquation(r)) + master_print(self.comm,"Reaction ", r+1, ": ", interface.getReactionEquation(r)) self.species_densities = np.array(interface.HostRead(libtps.t2bIndex.SpeciesDensities), copy=False) self.efield = np.array(interface.HostRead(libtps.t2bIndex.ElectricField), copy=False) self.heavy_temperature = np.array(interface.HostRead(libtps.t2bIndex.HeavyTemperature), copy=False) efieldAngularFreq = interface.EfieldAngularFreq() - print("Electric field angular frequency: ", efieldAngularFreq) + master_print(self.comm,"Electric field angular frequency: ", efieldAngularFreq) @@ -93,14 +99,14 @@ def push(self, interface): tps.chooseSolver() tps.initialize() -boltzmann = NullSolver() +boltzmann = NullSolver(comm) interface = libtps.Tps2Boltzmann(tps) tps.initInterface(interface) it = 0 max_iters = tps.getRequiredInput("cycle-avg-joule-coupled/max-iters") -print("Max Iters: ", max_iters) +master_print(comm,"Max Iters: ", max_iters) tps.solveBegin() while it < max_iters: @@ -113,7 +119,7 @@ def push(self, interface): tps.fetch(interface) it = it+1 - print("it, ", it) + master_print(comm, "it, ", it) tps.solveEnd() diff --git a/src/tps2Boltzmann.cpp b/src/tps2Boltzmann.cpp index 8de6dc79d..c646ee75e 100644 --- a/src/tps2Boltzmann.cpp +++ b/src/tps2Boltzmann.cpp @@ -118,7 +118,6 @@ int Tps2Boltzmann::_countBTEReactions() { } void Tps2Boltzmann::init(TPS::PlasmaSolver *flowSolver) { - std::cout << "Tps2Boltzmann::init is called" << std::endl; mfem::ParMesh *pmesh(flowSolver->getMesh()); if (use_h1fec_) fec_ = new mfem::H1_FECollection(order_, pmesh->Dimension(), basis_type_); From 0b57d65094cf6e8937b8f8a8ace4c8f56d369257 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Wed, 19 Feb 2025 15:37:58 -0600 Subject: [PATCH 08/32] Parallel bugfixes to set reaction rates --- src/reactingFlow.cpp | 17 +++++++++++------ src/reactingFlow.hpp | 3 ++- src/reaction.cpp | 3 ++- 3 files changed, 15 insertions(+), 8 deletions(-) diff --git a/src/reactingFlow.cpp b/src/reactingFlow.cpp index 015fcddc5..217246666 100644 --- a/src/reactingFlow.cpp +++ b/src/reactingFlow.cpp @@ -2887,16 +2887,21 @@ void ReactingFlow::push(TPS::Tps2Boltzmann &interface) { } void ReactingFlow::fetch(TPS::Tps2Boltzmann &interface) { + //NOTE: Chemistry is only address by truedofs indexes. mfem::ParFiniteElementSpace *reaction_rates_fes(&(interface.NativeFes(TPS::Tps2Boltzmann::Index::ReactionRates))); - externalReactionRates_.reset(new mfem::ParGridFunction(reaction_rates_fes)); - interface.interpolateToNativeFES(*externalReactionRates_, TPS::Tps2Boltzmann::Index::ReactionRates); + externalReactionRates_gf_.reset(new mfem::ParGridFunction(reaction_rates_fes)); + interface.interpolateToNativeFES(*externalReactionRates_gf_, TPS::Tps2Boltzmann::Index::ReactionRates); + externalReactionRates_gf_->GetTrueDofs(externalReactionRates_); + int size = sDofInt_; #if defined(_CUDA_) || defined(_HIP_) - const double *data(externalReactionRates_->Read()); - int size(externalReactionRates_->FESpace()->GetNDofs()); - assert(externalReactionRates_->FESpace()->GetOrdering() == mfem::Ordering::byNODES); + const double *data(externalReactionRates_.Read()); + //int size(externalReactionRates_->FESpace()->GetNDofs()); + assert(externalReactionRates_gf_->FESpace()->GetOrdering() == mfem::Ordering::byNODES); gpu::deviceSetChemistryReactionData<<<1, 1>>>(data, size, chemistry_); #else - chemistry_->setGridFunctionRates(*externalReactionRates_); + //chemistry_->setGridFunctionRates(*externalReactionRates_gf_); + const double *data(externalReactionRates_.HostRead()); + chemistry_->setRates(data, size); #endif } diff --git a/src/reactingFlow.hpp b/src/reactingFlow.hpp index a4219bb50..a70d6f8eb 100644 --- a/src/reactingFlow.hpp +++ b/src/reactingFlow.hpp @@ -104,7 +104,8 @@ class ReactingFlow : public ThermoChemModelBase { ArgonMixtureTransport *transport_ = NULL; Chemistry *chemistry_ = NULL; // External reaction rates when chemistry is implemented using the BTE option - std::unique_ptr externalReactionRates_; + std::unique_ptr externalReactionRates_gf_; //Has repeated interface dofs. + Vector externalReactionRates_; //Only true data std::vector speciesNames_; std::map atomMap_; diff --git a/src/reaction.cpp b/src/reaction.cpp index d119f58be..ef4db5d19 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -88,8 +88,8 @@ MFEM_HOST_DEVICE GridFunctionReaction::GridFunctionReaction(int comp) MFEM_HOST_DEVICE GridFunctionReaction::~GridFunctionReaction() {} MFEM_HOST_DEVICE void GridFunctionReaction::setData(const double *data, int size) { - data_ = data + comp_ * size_; size_ = size; + data_ = data + comp_ * size_; } void GridFunctionReaction::setGridFunction(const mfem::GridFunction &f) { @@ -103,6 +103,7 @@ void GridFunctionReaction::setGridFunction(const mfem::GridFunction &f) { #endif } + MFEM_HOST_DEVICE double GridFunctionReaction::computeRateCoefficient([[maybe_unused]] const double &T_h, [[maybe_unused]] const double &T_e, const int &dofindex, From b64923b0ab6f924154eb5f99d6fc4b81b79c1978 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Wed, 19 Feb 2025 16:18:46 -0600 Subject: [PATCH 09/32] Astyle --- src/reactingFlow.cpp | 2 +- src/reactingFlow.hpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/reactingFlow.cpp b/src/reactingFlow.cpp index 217246666..a544bf61c 100644 --- a/src/reactingFlow.cpp +++ b/src/reactingFlow.cpp @@ -2844,7 +2844,7 @@ double species_uniform(const Vector &coords, double t) { void ReactingFlow::push(TPS::Tps2Boltzmann &interface) { assert(interface.IsInitialized()); - //const int nscalardofs(vfes_->GetNDofs()); + // const int nscalardofs(vfes_->GetNDofs()); mfem::ParGridFunction *species = new mfem::ParGridFunction(&interface.NativeFes(TPS::Tps2Boltzmann::Index::SpeciesDensities)); diff --git a/src/reactingFlow.hpp b/src/reactingFlow.hpp index a70d6f8eb..a75746293 100644 --- a/src/reactingFlow.hpp +++ b/src/reactingFlow.hpp @@ -104,8 +104,8 @@ class ReactingFlow : public ThermoChemModelBase { ArgonMixtureTransport *transport_ = NULL; Chemistry *chemistry_ = NULL; // External reaction rates when chemistry is implemented using the BTE option - std::unique_ptr externalReactionRates_gf_; //Has repeated interface dofs. - Vector externalReactionRates_; //Only true data + std::unique_ptr externalReactionRates_gf_; // Has repeated interface dofs. + Vector externalReactionRates_; // Only true data std::vector speciesNames_; std::map atomMap_; From ef6d3550499273d7faed62bc21a613e1fc96b5d2 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Wed, 19 Feb 2025 16:25:47 -0600 Subject: [PATCH 10/32] Astyle --- src/reactingFlow.cpp | 6 +++--- src/reactingFlow.hpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/reactingFlow.cpp b/src/reactingFlow.cpp index a544bf61c..36a24a1d4 100644 --- a/src/reactingFlow.cpp +++ b/src/reactingFlow.cpp @@ -2887,7 +2887,7 @@ void ReactingFlow::push(TPS::Tps2Boltzmann &interface) { } void ReactingFlow::fetch(TPS::Tps2Boltzmann &interface) { - //NOTE: Chemistry is only address by truedofs indexes. + // NOTE: Chemistry is only address by truedofs indexes. mfem::ParFiniteElementSpace *reaction_rates_fes(&(interface.NativeFes(TPS::Tps2Boltzmann::Index::ReactionRates))); externalReactionRates_gf_.reset(new mfem::ParGridFunction(reaction_rates_fes)); interface.interpolateToNativeFES(*externalReactionRates_gf_, TPS::Tps2Boltzmann::Index::ReactionRates); @@ -2895,11 +2895,11 @@ void ReactingFlow::fetch(TPS::Tps2Boltzmann &interface) { int size = sDofInt_; #if defined(_CUDA_) || defined(_HIP_) const double *data(externalReactionRates_.Read()); - //int size(externalReactionRates_->FESpace()->GetNDofs()); + // int size(externalReactionRates_->FESpace()->GetNDofs()); assert(externalReactionRates_gf_->FESpace()->GetOrdering() == mfem::Ordering::byNODES); gpu::deviceSetChemistryReactionData<<<1, 1>>>(data, size, chemistry_); #else - //chemistry_->setGridFunctionRates(*externalReactionRates_gf_); + // chemistry_->setGridFunctionRates(*externalReactionRates_gf_); const double *data(externalReactionRates_.HostRead()); chemistry_->setRates(data, size); #endif diff --git a/src/reactingFlow.hpp b/src/reactingFlow.hpp index a75746293..56fc1bc35 100644 --- a/src/reactingFlow.hpp +++ b/src/reactingFlow.hpp @@ -104,8 +104,8 @@ class ReactingFlow : public ThermoChemModelBase { ArgonMixtureTransport *transport_ = NULL; Chemistry *chemistry_ = NULL; // External reaction rates when chemistry is implemented using the BTE option - std::unique_ptr externalReactionRates_gf_; // Has repeated interface dofs. - Vector externalReactionRates_; // Only true data + std::unique_ptr externalReactionRates_gf_; // Has repeated interface dofs. + Vector externalReactionRates_; // Only true data std::vector speciesNames_; std::map atomMap_; From 62550381828e36f006305c4dc190d959dab227c8 Mon Sep 17 00:00:00 2001 From: milindasf Date: Wed, 26 Feb 2025 23:34:23 -0600 Subject: [PATCH 11/32] lomach bte 4sp --- src/tps-bte_0d3v.py | 45 ++++++++++++++++++++++++--------------------- 1 file changed, 24 insertions(+), 21 deletions(-) diff --git a/src/tps-bte_0d3v.py b/src/tps-bte_0d3v.py index bc1976955..ce7d65bd5 100755 --- a/src/tps-bte_0d3v.py +++ b/src/tps-bte_0d3v.py @@ -69,7 +69,8 @@ def min_mean_max(a, comm: MPI.Comm): from parla import Parla from parla.tasks import spawn, TaskSpace from parla.devices import cpu, gpu - except: + except Exception as e: + print(e) print("Error occured during Parla import. Please make sure Parla is installed properly.") sys.exit(0) @@ -122,7 +123,7 @@ class BoltzmannSolverParams(): n0 = 3.22e22 #[m^{-3}] rand_seed = 0 - use_clstr_inp = False + use_clstr_inp = True clstr_maxiter = 10 clstr_threshold = 1e-3 @@ -132,27 +133,20 @@ class TPSINDEX(): """ simple index map to differnt fields, from the TPS arrays """ - # ION_IDX = 0 # ion density index - # ELE_IDX = 1 # electron density index - # NEU_IDX = 2 # neutral density index - EF_RE_IDX = 0 # Re(E) index EF_IM_IDX = 1 # Im(E) index # in future we need to setup this methodically # here key denotes the idx running from 0, nreactions-1 # value denotes the reaction index in the qoi array - RR_IDX = {0 : 4 , 1 : 5 , 2 : 6, 3 : 7, 4 : 1 , 5 : 2, 6 : 3 } - + RR_IDX = {0 : 2 , 1 : 3 , 2 : 1} - ION_IDX = 3 - ELE_IDX = 4 - NEU_IDX = 5 + ION_IDX = 1 + ELE_IDX = 2 + NEU_IDX = 3 EX1_IDX = 0 - EX2_IDX = 1 - EX3_IDX = 2 - MOLE_FRAC_IDX = {0: NEU_IDX, 1: EX1_IDX , 2: EX2_IDX , 3: EX3_IDX} + MOLE_FRAC_IDX = {0: NEU_IDX, 1: EX1_IDX} def k_means(x, num_clusters, xi=None, max_iter=1000, thresh=1e-12, rand_seed=0, xp=np): assert x.ndim == 2, "observations must me 2d array" @@ -344,7 +338,7 @@ def grid_setup(self, interface): # active_grid_idx.append(grid_idx) # self.active_grid_idx = active_grid_idx #[i for i in range(self.param.n_grids)] - self.active_grid_idx = [2,3]#[i for i in range(self.param.n_grids)] + self.active_grid_idx = [i for i in range(self.param.n_grids)] self.sub_clusters_run = False return @@ -376,9 +370,10 @@ def solve_wo_parla(self): ff , qoi = self.bte_solver.solve(grid_idx, f0, self.param.atol, self.param.rtol, self.param.max_iter, self.param.solver_type) self.qoi[grid_idx] = qoi self.ff [grid_idx] = ff - except: + except Exception as e: + print(e) print("solver failed for v-space gird no %d"%(grid_idx), flush=True) - sys.exit(-1) + self.comm.Abort(0) else: with xp.cuda.Device(dev_id): try: @@ -386,9 +381,10 @@ def solve_wo_parla(self): ff , qoi = self.bte_solver.solve(grid_idx, f0, self.param.atol, self.param.rtol, self.param.max_iter, self.param.solver_type) self.qoi[grid_idx] = qoi self.ff [grid_idx] = ff - except: + except Exception as e: + print(e) print("solver failed for v-space gird no %d"%(grid_idx), flush=True) - sys.exit(-1) + self.comm.Abort(0) t2 = time() print("time for boltzmann v-space solve = %.4E"%(t2- t1), flush=True) @@ -830,15 +826,22 @@ def solve(self): dev_id = gidx_to_device_map(grid_idx,n_grids) def t1(): + # print("rank [%d/%d] BTE launching grid %d on %s"%(rank, npes, grid_idx, dev_id), flush=True) + # f0 = self.bte_solver.get_boltzmann_parameter(grid_idx, "u0") + # ff , qoi = self.bte_solver.solve(grid_idx, f0, self.param.atol, self.param.rtol, self.param.max_iter, self.param.solver_type) + # self.ff[grid_idx] = ff + # self.qoi[grid_idx] = qoi try: print("rank [%d/%d] BTE launching grid %d on %s"%(rank, npes, grid_idx, dev_id), flush=True) f0 = self.bte_solver.get_boltzmann_parameter(grid_idx, "u0") ff , qoi = self.bte_solver.solve(grid_idx, f0, self.param.atol, self.param.rtol, self.param.max_iter, self.param.solver_type) self.ff[grid_idx] = ff self.qoi[grid_idx] = qoi - except: + except Exception as e: + print(e) print("rank [%d/%d] solver failed for v-space gird no %d"%(self.rankG, self.npesG, grid_idx), flush=True) - sys.exit(-1) + self.comm.Abort(0) + #sys.exit(-1) with cp.cuda.Device(dev_id): t1() From a39871b99b4781a5038a3329e2ede61b3c78cc52 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Fri, 28 Feb 2025 10:00:37 -0600 Subject: [PATCH 12/32] Modify the mock solver to handle more than one reaction --- src/tps-time-loop.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/tps-time-loop.py b/src/tps-time-loop.py index cd5018298..0ae567f4d 100755 --- a/src/tps-time-loop.py +++ b/src/tps-time-loop.py @@ -40,7 +40,10 @@ def solve(self): def push(self, interface): rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False) - rates[:] = 0. + + n_reactions =interface.nComponents(libtps.t2bIndex.ReactionRates) + for r in range(n_reactions): + rates[r*self.heavy_temperature.shape[0]:(r+1)*self.heavy_temperature.shape[0]] = (10.**(-r))*1e-6*self.heavy_temperature class ArrheniusSolver: def __init__(self,comm): From cdd6e561bda017360add7b64263c7285059b3c50 Mon Sep 17 00:00:00 2001 From: milindasf Date: Fri, 21 Mar 2025 17:28:29 -0500 Subject: [PATCH 13/32] 4sp with deexcitation --- src/tps-bte_0d3v.py | 36 +++++++++++++++++++++++++++--------- src/tps.py | 11 ++++++++++- 2 files changed, 37 insertions(+), 10 deletions(-) diff --git a/src/tps-bte_0d3v.py b/src/tps-bte_0d3v.py index ce7d65bd5..594cff330 100755 --- a/src/tps-bte_0d3v.py +++ b/src/tps-bte_0d3v.py @@ -15,6 +15,7 @@ import scipy.cluster import threading import datetime +import h5py # Use asynchronous stream ordered memory #cp.cuda.set_allocator(cp.cuda.MemoryAsyncPool().malloc) @@ -139,7 +140,7 @@ class TPSINDEX(): # in future we need to setup this methodically # here key denotes the idx running from 0, nreactions-1 # value denotes the reaction index in the qoi array - RR_IDX = {0 : 2 , 1 : 3 , 2 : 1} + RR_IDX = {0 : 2 , 1 : 4 , 2 : 1, 3 : 3} ION_IDX = 1 ELE_IDX = 2 @@ -858,7 +859,7 @@ def push(self, interface): tps_npts = len(heavy_temp) n_reactions = interface.nComponents(libtps.t2bIndex.ReactionRates) - rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False).reshape((n_reactions, tps_npts)) + rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False).reshape((n_reactions, tps_npts), copy=False) if (use_interp==True): if(n_reactions>0): @@ -886,9 +887,17 @@ def t1(): with cp.cuda.Device(dev_id): t1() - - rates = rates.reshape((-1)) + + rates = rates.reshape((-1), copy=False) rates[rates<0] = 0.0 + + # fname = self.param.out_fname+"_push_rank_%d_npes_%d.h5"%(self.rankG, self.npesG) + # with h5py.File(fname, 'w') as F: + # F.create_dataset("Tg[K]" , data=heavy_temp) + # F.create_dataset("rates" , data=rates) + # F.create_dataset("rates2" , data=np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False)) + # F.close() + else: if(n_reactions>0): rates[:,:] = 0.0 @@ -909,7 +918,7 @@ def t1(): with cp.cuda.Device(dev_id): t1() - rates = rates.reshape((-1)) + rates = rates.reshape((-1), copy=False) rates[rates<0] = 0.0 return @@ -1284,7 +1293,7 @@ async def push_async(self, interface): tps_npts = len(heavy_temp) n_reactions = interface.nComponents(libtps.t2bIndex.ReactionRates) - rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False).reshape((n_reactions, tps_npts)) + rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False).reshape((n_reactions, tps_npts), copy=False) if (use_interp==True): if(n_reactions>0): @@ -1819,8 +1828,8 @@ def __main__(): print("tps steps per cycle : ", tps_sper_cycle, "bte_steps per cycle", bte_sper_cycle, flush=True) while (iter Date: Wed, 26 Mar 2025 00:52:57 -0500 Subject: [PATCH 14/32] tempory remove tps(comm) and additional copy=False in reshape function to be executed on frontera --- src/tps-bte_0d3v.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/tps-bte_0d3v.py b/src/tps-bte_0d3v.py index 594cff330..e5efd6b9b 100755 --- a/src/tps-bte_0d3v.py +++ b/src/tps-bte_0d3v.py @@ -859,7 +859,7 @@ def push(self, interface): tps_npts = len(heavy_temp) n_reactions = interface.nComponents(libtps.t2bIndex.ReactionRates) - rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False).reshape((n_reactions, tps_npts), copy=False) + rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False).reshape((n_reactions, tps_npts)) if (use_interp==True): if(n_reactions>0): @@ -888,7 +888,7 @@ def t1(): with cp.cuda.Device(dev_id): t1() - rates = rates.reshape((-1), copy=False) + rates = rates.reshape((-1)) rates[rates<0] = 0.0 # fname = self.param.out_fname+"_push_rank_%d_npes_%d.h5"%(self.rankG, self.npesG) @@ -918,7 +918,7 @@ def t1(): with cp.cuda.Device(dev_id): t1() - rates = rates.reshape((-1), copy=False) + rates = rates.reshape((-1)) rates[rates<0] = 0.0 return @@ -1293,7 +1293,7 @@ async def push_async(self, interface): tps_npts = len(heavy_temp) n_reactions = interface.nComponents(libtps.t2bIndex.ReactionRates) - rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False).reshape((n_reactions, tps_npts), copy=False) + rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False).reshape((n_reactions, tps_npts)) if (use_interp==True): if(n_reactions>0): @@ -1781,7 +1781,8 @@ def driver_wo_parla(comm): def __main__(): # TPS solver profile_tt[pp.TPS_SETUP].start() - tps = libtps.Tps(comm) + #tps = libtps.Tps(comm) + tps = libtps.Tps() tps.parseCommandLineArgs(sys.argv) tps.parseInput() tps.chooseDevices() From 8105b710c24b48b9e55050b545eb9903ef39f992 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Fri, 2 May 2025 21:05:11 -0500 Subject: [PATCH 15/32] Exchange tabulated rates from python -- hard-coded --- src/tps-time-loop.py | 55 ++++++++++++++++++++++++++------------------ 1 file changed, 33 insertions(+), 22 deletions(-) diff --git a/src/tps-time-loop.py b/src/tps-time-loop.py index 0ae567f4d..aa6e40898 100755 --- a/src/tps-time-loop.py +++ b/src/tps-time-loop.py @@ -2,6 +2,9 @@ import sys import os import numpy as np +import h5py +import scipy +import scipy.interpolate from mpi4py import MPI @@ -45,20 +48,31 @@ def push(self, interface): for r in range(n_reactions): rates[r*self.heavy_temperature.shape[0]:(r+1)*self.heavy_temperature.shape[0]] = (10.**(-r))*1e-6*self.heavy_temperature -class ArrheniusSolver: - def __init__(self,comm): + +class TabulatedSolver: + def __init__(self, comm): self.comm = comm - self.UNIVERSALGASCONSTANT = 8.3144598; # J * mol^(-1) * K^(-1) self.species_densities = None self.efield = None self.heavy_temperature = None - self.reaction_rates = [None, None] - #Reaction 1: 'Ar + E => Ar.+1 + 2 E', - #Reaction 2: 'Ar.+1 + 2 E => Ar + E' - self.A = [74072.331348, 5.66683445516e-20] - self.b = [1.511, 0.368] - self.E = [1176329.772504, -377725.908714] # [J/mol] + self.tables = self._read_tables() + self.rates = [] + + def _read_tables(self): + filenames = ["./rate-coefficients/Ionization_Ground.h5", + "./rate-coefficients/Ionization_Lumped.h5", + "./rate-coefficients/Excitation_Lumped.h5"] + tables = [] + for filename in filenames: + with h5py.File(filename, 'r') as fid: + Tcoeff = fid['table'][:] + + tables.append(scipy.interpolate.interp1d(Tcoeff[:,0], Tcoeff[:,1], kind='linear', + bounds_error=False, fill_value='extrapolate')) + + return tables + def fetch(self, interface): n_reactions =interface.nComponents(libtps.t2bIndex.ReactionRates) for r in range(n_reactions): @@ -70,20 +84,17 @@ def fetch(self, interface): efieldAngularFreq = interface.EfieldAngularFreq() master_print(self.comm,"Electric field angular frequency: ", efieldAngularFreq) - - def solve(self): - #A_ * pow(temp, b_) * exp(-E_ / UNIVERSALGASCONSTANT / temp); - self.reaction_rates = [A * np.power(self.heavy_temperature, b) * - np.exp(-E/(self.UNIVERSALGASCONSTANT * self.heavy_temperature)) - for A,b,E in zip(self.A, self.b, self.E) ] + self.rates = [] + for table in self.tables: + self.rates.append(table(self.heavy_temperature)) def push(self, interface): - n_reactions =interface.nComponents(libtps.t2bIndex.ReactionRates) - if n_reactions >= 2: - rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False) - rates[0:self.heavy_temperature.shape[0]] = self.reaction_rates[0] - rates[self.heavy_temperature.shape[0]:] = self.reaction_rates[1] + rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False) + offset = 0 + for rate in self.rates: + rates[offset:offset+rate.shape[0]] = rate + offset = offset+rate.shape[0] @@ -102,7 +113,7 @@ def push(self, interface): tps.chooseSolver() tps.initialize() -boltzmann = NullSolver(comm) +boltzmann = TabulatedSolver(comm) interface = libtps.Tps2Boltzmann(tps) tps.initInterface(interface) @@ -113,13 +124,13 @@ def push(self, interface): tps.solveBegin() while it < max_iters: - tps.solveStep() tps.push(interface) boltzmann.fetch(interface) boltzmann.solve() boltzmann.push(interface) interface.saveDataCollection(cycle=it, time=it) tps.fetch(interface) + tps.solveStep() it = it+1 master_print(comm, "it, ", it) From 7a9922351f2160ea6d6d022418eeabfa1c463e32 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Tue, 6 May 2025 16:57:21 -0500 Subject: [PATCH 16/32] Let tps-time-loop parse input file --- src/tps-time-loop.py | 39 ++++++++++++++++++++++++++++++++++----- 1 file changed, 34 insertions(+), 5 deletions(-) diff --git a/src/tps-time-loop.py b/src/tps-time-loop.py index aa6e40898..1b9f6a40a 100755 --- a/src/tps-time-loop.py +++ b/src/tps-time-loop.py @@ -6,6 +6,8 @@ import scipy import scipy.interpolate +import configparser + from mpi4py import MPI def master_print(comm: MPI.Comm, *args, **kwargs) -> None: @@ -50,8 +52,9 @@ def push(self, interface): class TabulatedSolver: - def __init__(self, comm): + def __init__(self, comm, config): self.comm = comm + self.config = config self.species_densities = None self.efield = None self.heavy_temperature = None @@ -59,10 +62,23 @@ def __init__(self, comm): self.tables = self._read_tables() self.rates = [] + def _findPythonReactions(self): + filenames = [] + nreactions = self.config.getint("reactions","number_of_reactions",fallback=0) + for ir in range(nreactions): + sublist = self.config["reactions/reaction{0:d}".format(ir+1)] + rtype = sublist["model"] + if rtype == "bte": + filenames.append(sublist["tabulated/filename"].strip("'")) + + return filenames + def _read_tables(self): - filenames = ["./rate-coefficients/Ionization_Ground.h5", - "./rate-coefficients/Ionization_Lumped.h5", - "./rate-coefficients/Excitation_Lumped.h5"] + filenames = self._findPythonReactions() + + #["./rate-coefficients/Ionization_Ground.h5", + # "./rate-coefficients/Ionization_Lumped.h5", + # "./rate-coefficients/Excitation_Lumped.h5"] tables = [] for filename in filenames: with h5py.File(filename, 'r') as fid: @@ -107,13 +123,26 @@ def push(self, interface): # TPS solver tps = libtps.Tps(comm) + + tps.parseCommandLineArgs(sys.argv) tps.parseInput() tps.chooseDevices() tps.chooseSolver() tps.initialize() -boltzmann = TabulatedSolver(comm) +ini_name = '' +if '-run' in sys.argv: + ini_name = sys.argv[sys.argv.index('-run') + 1 ] +elif '-runFile' in sys.argv: + ini_name = sys.argv[sys.argv.index('-runFile') + 1 ] +else: + exit(-1) + +config = configparser.ConfigParser() +config.read(ini_name) + +boltzmann = TabulatedSolver(comm, config) interface = libtps.Tps2Boltzmann(tps) tps.initInterface(interface) From d38e03a63f441306124e9c2f5a905c032809efc7 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Tue, 6 May 2025 17:09:30 -0500 Subject: [PATCH 17/32] Bug fix --- src/tps-time-loop.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/tps-time-loop.py b/src/tps-time-loop.py index 1b9f6a40a..a37c48156 100755 --- a/src/tps-time-loop.py +++ b/src/tps-time-loop.py @@ -134,11 +134,13 @@ def push(self, interface): ini_name = '' if '-run' in sys.argv: ini_name = sys.argv[sys.argv.index('-run') + 1 ] -elif '-runFile' in sys.argv: - ini_name = sys.argv[sys.argv.index('-runFile') + 1 ] +elif '--runFile' in sys.argv: + ini_name = sys.argv[sys.argv.index('--runFile') + 1 ] else: + print("Could not parse command line in python. GOOD BYE!") exit(-1) +print(ini_name) config = configparser.ConfigParser() config.read(ini_name) From e860e30e9774e0db65a6c783a7a8e8fc14fb830c Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Wed, 14 May 2025 11:24:37 -0500 Subject: [PATCH 18/32] Fix tps.py test --- src/tps.py | 11 ++--------- test/cyl3d.python.test | 35 ++++++++++++++++++++++++++++++++--- 2 files changed, 34 insertions(+), 12 deletions(-) diff --git a/src/tps.py b/src/tps.py index 0c7f817bc..780bb8c79 100755 --- a/src/tps.py +++ b/src/tps.py @@ -18,15 +18,8 @@ tps.chooseDevices() tps.chooseSolver() tps.initialize() -tps.solveStep() -#tps.solve() +tps.solve() + -for i in range(20): - comm.Barrier() - t1 = time() - tps.solveStep() - comm.Barrier() - t2 = time() - print("step [%04d] solve time (s)= %.4E"%(i, t2-t1)) sys.exit (tps.getStatus()) diff --git a/test/cyl3d.python.test b/test/cyl3d.python.test index 8036c1a96..23be6c507 100755 --- a/test/cyl3d.python.test +++ b/test/cyl3d.python.test @@ -9,6 +9,8 @@ EXE="../src/tps.py" setup() { SOLN_FILE=restart_output.sol.h5 REF_FILE=ref_solns/cyl3d_coarse.4iters.h5 + rm -rf planeData + mkdir planeData } @test "[$TEST] check for input file $RUNFILE" { @@ -109,7 +111,7 @@ setup() { ./soln_differ $SOLN_FILE $REF_FILE } -@test "[$TEST] verify serial restart consistency using 2 mpi tasks" { +@test "[$TEST] verify singleFileWrite restart consistency using 2 mpi tasks" { RUNFILE=inputs/input.4iters.cyl.ini RUNFILE_PART=${RUNFILE}.part SOLN_FILE=restart_output.sol.h5 @@ -136,6 +138,33 @@ setup() { rm -f restart_output.sol.*.h5 } +@test "[$TEST] verify singleFileReadWRite restart consistency using 2 mpi tasks" { + RUNFILE=inputs/input.4iters.cyl.ini + RUNFILE_PART=${RUNFILE}.part + SOLN_FILE=restart_output.sol.h5 + + rm -f $RUNFILE_PART + rm -f $SOLN_FILE + + # copy the run file, b/c I want to edit it + test -s $RUNFILE + cp $RUNFILE $RUNFILE_PART + + # request single file restart + echo "[io]" >> $RUNFILE_PART + echo "restartMode = singleFileReadWrite" >> $RUNFILE_PART + + # run 2 tasks + mpirun -np 2 $EXE --runFile $RUNFILE_PART + test -s $SOLN_FILE_SERIAL + + ./soln_differ $SOLN_FILE $REF_FILE + + # clean up if we got this far + rm -f $RUNFILE_PART + rm -f restart_output.sol.*.h5 +} + @test "[$TEST] verify serial restart consistency using 4 mpi tasks" { RUNFILE=inputs/input.4iters.cyl.ini RUNFILE_PART=${RUNFILE}.part @@ -194,7 +223,7 @@ setup() { cp $RUNFILE $RUNFILE_MOD echo "[io]" >> $RUNFILE_MOD echo "restartMode = singleFileRead" >> $RUNFILE_MOD - echo "enableRestart = true" >> $RUNFILE_MOD + echo "enableRestart = True" >> $RUNFILE_MOD mpirun -np 4 $EXE --runFile $RUNFILE_MOD test -s partition.4p.h5 @@ -256,7 +285,7 @@ setup() { cp $RUNFILE $RUNFILE_MOD echo "[io]" >> $RUNFILE_MOD echo "restartMode = singleFileWrite" >> $RUNFILE_MOD - echo "enableRestart = true" >> $RUNFILE_MOD + echo "enableRestart = True" >> $RUNFILE_MOD mpirun -np 2 $EXE --runFile $RUNFILE_MOD test -s partition.2p.h5 From 546f206b9a5747f25e1592b5e90be7515c0d605f Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Wed, 14 May 2025 12:55:40 -0500 Subject: [PATCH 19/32] Update Dockerfile Adding scipy to the contaneir --- docker/test/Dockerfile | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docker/test/Dockerfile b/docker/test/Dockerfile index d3869b629..af28d5b97 100644 --- a/docker/test/Dockerfile +++ b/docker/test/Dockerfile @@ -147,11 +147,13 @@ RUN pip3 install gcovr RUN yum -y install gsl-gnu9-ohpc -RUN pip3 install "pybind11[global]" - RUN . /etc/profile.d/lmod.sh \ + && pip3 install "numpy<2.0" \ + && pip3 install scipy \ && pip3 install mpi4py +RUN pip3 install "pybind11[global]" + # Register new libs installed into /usr/local/lib with linker RUN echo "/usr/local/lib" > /etc/ld.so.conf.d/class.conf RUN ldconfig From 7ce879c85df42d6af6412ad94b46370a2a45099e Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Thu, 15 May 2025 13:25:01 -0500 Subject: [PATCH 20/32] Some docker fixes --- docker/test/Dockerfile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docker/test/Dockerfile b/docker/test/Dockerfile index af28d5b97..6f64c5a8a 100644 --- a/docker/test/Dockerfile +++ b/docker/test/Dockerfile @@ -29,7 +29,6 @@ RUN yum -y install valgrind-ohpc RUN yum -y install boost-gnu9-mpich-ohpc RUN yum -y install python38-pip -RUN pip3 install matplotlib RUN yum -y install procps-ng RUN yum -y install diffutils @@ -148,8 +147,9 @@ RUN pip3 install gcovr RUN yum -y install gsl-gnu9-ohpc RUN . /etc/profile.d/lmod.sh \ - && pip3 install "numpy<2.0" \ + && pip3 install "numpy<2" \ && pip3 install scipy \ + && pip3 install matplotlib \ && pip3 install mpi4py RUN pip3 install "pybind11[global]" From 5c0a6772aa92fada6a03486668931306dc6cbac2 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Thu, 15 May 2025 14:50:23 -0500 Subject: [PATCH 21/32] Iron out the Docker container --- docker/test/Dockerfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docker/test/Dockerfile b/docker/test/Dockerfile index 6f64c5a8a..f8c95e990 100644 --- a/docker/test/Dockerfile +++ b/docker/test/Dockerfile @@ -150,7 +150,7 @@ RUN . /etc/profile.d/lmod.sh \ && pip3 install "numpy<2" \ && pip3 install scipy \ && pip3 install matplotlib \ - && pip3 install mpi4py + && pip3 install "mpi4py<4" RUN pip3 install "pybind11[global]" From fe0c2bdf13eb28aafa0b809ce298ea41a0e83c71 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Thu, 15 May 2025 16:02:54 -0500 Subject: [PATCH 22/32] Add h5py to the Docker container --- docker/test/Dockerfile | 1 + 1 file changed, 1 insertion(+) diff --git a/docker/test/Dockerfile b/docker/test/Dockerfile index f8c95e990..d22131292 100644 --- a/docker/test/Dockerfile +++ b/docker/test/Dockerfile @@ -150,6 +150,7 @@ RUN . /etc/profile.d/lmod.sh \ && pip3 install "numpy<2" \ && pip3 install scipy \ && pip3 install matplotlib \ + && pip3 install h5py \ && pip3 install "mpi4py<4" RUN pip3 install "pybind11[global]" From 04f0e834a2dad91fb40f42f65e71ce8796892888 Mon Sep 17 00:00:00 2001 From: milindasf Date: Fri, 16 May 2025 15:50:12 -0500 Subject: [PATCH 23/32] temperature based crs initial attempt --- src/tps-bte_0d3v.py | 113 +++++++++++++++++++++++--------------------- 1 file changed, 59 insertions(+), 54 deletions(-) diff --git a/src/tps-bte_0d3v.py b/src/tps-bte_0d3v.py index e5efd6b9b..91836a5bf 100755 --- a/src/tps-bte_0d3v.py +++ b/src/tps-bte_0d3v.py @@ -59,10 +59,18 @@ def min_mean_max(a, comm: MPI.Comm): # set path to C++ TPS library path = os.path.abspath(os.path.dirname(sys.argv[0])) sys.path.append(path + "/.libs") +sys.path.append(path + "/../../torch-chemistry") +import argon.scripts.synthetic_cs as synthetic_cs sys.path.append(path + "/../../boltzmann/BESolver/python") + +crs_folder = path + "/../../torch-chemistry/argon/scripts/crs_lxcat" +if not os.path.exists(crs_folder): + os.makedirs(crs_folder) + import libtps from bte_0d3v_batched import bte_0d3v_batched as BoltzmannSolver import utils as bte_utils +print(bte_utils) WITH_PARLA = 0 if WITH_PARLA: @@ -134,6 +142,10 @@ class TPSINDEX(): """ simple index map to differnt fields, from the TPS arrays """ + # ION_IDX = 0 # ion density index + # ELE_IDX = 1 # electron density index + # NEU_IDX = 2 # neutral density index + EF_RE_IDX = 0 # Re(E) index EF_IM_IDX = 1 # Im(E) index @@ -314,7 +326,15 @@ def grid_setup(self, interface): Te = xp.array([Te_b[b_idx] for b_idx in range(self.param.n_grids)]) # xp.ones(self.param.n_grids) * self.param.Te vth = np.sqrt(2* self.param.kB * Te * self.param.ev_to_K /self.param.me) ev_max = (6 * vth / self.param.c_gamma)**2 - self.bte_solver = BoltzmannSolver(self.param, ev_max , Te , nr, lm_modes, self.param.n_grids, self.param.collisions) + + # generate crs Tg depended crs data + col_cs = list() + for idx in range(self.param.n_grids): + cs_fname = "%s/crs_rank_%06d_npes_%06d_%06d.txt"%(crs_folder,self.rankG, self.npesG, idx) + synthetic_cs.gen_lxcat_file(cs_fname, Te_b[idx] * self.param.ev_to_K) + col_cs.append(cs_fname) + + self.bte_solver = BoltzmannSolver(self.param, ev_max , Te , nr, lm_modes, self.param.n_grids, col_cs) # compute BTE operators for grid_idx in range(self.param.n_grids): @@ -354,13 +374,15 @@ def solve_wo_parla(self): n_grids = self.param.n_grids gidx_to_device_map = self.gidx_to_device_map - self.qoi = [None for grid_idx in range(self.param.n_grids)] - self.ff = [None for grid_idx in range(self.param.n_grids)] - coll_list = self.bte_solver.get_collision_list() - coll_names = self.bte_solver.get_collision_names() + self.qoi = [None for grid_idx in range(self.param.n_grids)] + self.ff = [None for grid_idx in range(self.param.n_grids)] + if csv_write: - data_csv = np.empty((self.tps_npts, 8 + len(coll_list))) + for grid_idx in range(n_grids): + assert len(self.bte_solver.get_collision_list()[i]) == len(self.bte_solver.get_collision_list()[0]) + + data_csv = np.empty((self.tps_npts, 8 + len(self.bte_solver.get_collision_list()[0]))) t1 = time() for grid_idx in range(n_grids): @@ -392,6 +414,10 @@ def solve_wo_parla(self): if (self.param.export_csv ==1 or self.param.plot_data==1): for grid_idx in range(n_grids): + + coll_list = self.bte_solver.get_collision_list()[grid_idx] + coll_names = self.bte_solver.get_collision_names()[grid_idx] + dev_id = gidx_to_device_map(grid_idx, n_grids) ff = self.ff[grid_idx] qoi = self.qoi[grid_idx] @@ -480,7 +506,9 @@ def fetch(self, interface): efield = np.array(interface.HostRead(libtps.t2bIndex.ElectricField), copy=False).reshape((2, tps_npts)) species_densities = np.array(interface.HostRead(libtps.t2bIndex.SpeciesDensities), copy=False).reshape(nspecies, tps_npts) - cs_avail_species = self.bte_solver._avail_species + cs_avail_species = self.bte_solver._avail_species[0] + for i in range(len(self.bte_solver._avail_species)): + assert cs_avail_species == self.bte_solver._avail_species[i] n0 = np.sum(species_densities, axis=0) - species_densities[TPSINDEX.ELE_IDX] ns_by_n0 = np.concatenate([species_densities[TPSINDEX.MOLE_FRAC_IDX[i]]/n0 for i in range(len(cs_avail_species))]).reshape((len(cs_avail_species), tps_npts)) @@ -494,21 +522,23 @@ def fetch(self, interface): Ey = efield[1] ne = species_densities[TPSINDEX.ELE_IDX] - coll_list = self.bte_solver.get_collision_list() - coll_names = self.bte_solver.get_collision_names() - cs_data = self.bte_solver.get_cross_section_data() - cs_species = list() - for col_idx, (k,v) in enumerate(cs_data.items()): - cs_species.append(v["species"]) - - cs_species = list(sorted(set(cs_species), key=cs_species.index)) data_csv = np.concatenate([(Ex).reshape((-1, 1)), (Ey).reshape((-1, 1)), (Tg).reshape((-1, 1)), (ne/n0).reshape((-1, 1))] + [ns_by_n0[i].reshape((-1, 1)) for i in range(ns_by_n0.shape[0])] + [n0.reshape(-1, 1)], axis=1) for grid_idx in self.active_grid_idx: + coll_list = self.bte_solver.get_collision_list()[grid_idx] + coll_names = self.bte_solver.get_collision_names()[grid_idx] + cs_data = self.bte_solver.get_cross_section_data()[grid_idx] + + cs_species = list() + for col_idx, (k,v) in enumerate(cs_data.items()): + cs_species.append(v["species"]) + cs_species = list(sorted(set(cs_species), key=cs_species.index)) + + with open("%s/%s.csv"%(self.param.output_dir, "tps_fetch_grid_%02d_rank_%02d_npes_%02d"%(grid_idx, self.rankG, self.npesG)), 'w', encoding='UTF8') as f: writer = csv.writer(f,delimiter=',') # write the header @@ -817,21 +847,12 @@ def solve(self): self.qoi = [None for grid_idx in range(self.param.n_grids)] self.ff = [None for grid_idx in range(self.param.n_grids)] - coll_list = self.bte_solver.get_collision_list() - coll_names = self.bte_solver.get_collision_names() - - if csv_write: - data_csv = np.empty((self.tps_npts, 8 + len(coll_list))) - + + for grid_idx in self.active_grid_idx: dev_id = gidx_to_device_map(grid_idx,n_grids) def t1(): - # print("rank [%d/%d] BTE launching grid %d on %s"%(rank, npes, grid_idx, dev_id), flush=True) - # f0 = self.bte_solver.get_boltzmann_parameter(grid_idx, "u0") - # ff , qoi = self.bte_solver.solve(grid_idx, f0, self.param.atol, self.param.rtol, self.param.max_iter, self.param.solver_type) - # self.ff[grid_idx] = ff - # self.qoi[grid_idx] = qoi try: print("rank [%d/%d] BTE launching grid %d on %s"%(rank, npes, grid_idx, dev_id), flush=True) f0 = self.bte_solver.get_boltzmann_parameter(grid_idx, "u0") @@ -874,7 +895,6 @@ def t1(): qoi = self.bte_solver.compute_QoIs(grid_idx, h_curr, effective_mobility=False) rr_cpu = xp.asnumpy(qoi["rates"]) inp_mask = xp.asnumpy(self.sub_cluster_c_lbl[grid_idx]) == np.arange(self.param.n_sub_clusters)[:, None] - rr_interp = np.zeros((n_reactions, len(gidx_to_pidx_map[grid_idx]))) for c_idx in range(self.param.n_sub_clusters): @@ -890,14 +910,6 @@ def t1(): rates = rates.reshape((-1)) rates[rates<0] = 0.0 - - # fname = self.param.out_fname+"_push_rank_%d_npes_%d.h5"%(self.rankG, self.npesG) - # with h5py.File(fname, 'w') as F: - # F.create_dataset("Tg[K]" , data=heavy_temp) - # F.create_dataset("rates" , data=rates) - # F.create_dataset("rates2" , data=np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False)) - # F.close() - else: if(n_reactions>0): rates[:,:] = 0.0 @@ -934,7 +946,9 @@ async def fetch_asnyc(self, interface): efield = np.array(interface.HostRead(libtps.t2bIndex.ElectricField), copy=False).reshape((2, tps_npts)) species_densities = np.array(interface.HostRead(libtps.t2bIndex.SpeciesDensities), copy=False).reshape(nspecies, tps_npts) - cs_avail_species = self.bte_solver._avail_species + cs_avail_species = self.bte_solver._avail_species[0] + for i in range(len(self.bte_solver._avail_species)): + assert cs_avail_species == self.bte_solver._avail_species[i] n0 = np.sum(species_densities, axis=0) - species_densities[TPSINDEX.ELE_IDX] ns_by_n0 = np.concatenate([species_densities[TPSINDEX.MOLE_FRAC_IDX[i]]/n0 for i in range(len(cs_avail_species))]).reshape((len(cs_avail_species), tps_npts)) @@ -1252,16 +1266,13 @@ async def solve_async(self): self.ff = [None for grid_idx in range(self.param.n_grids)] num_gpus = len(gpu) assert num_gpus == self.num_gpus_per_node, "CuPy and Parla number of GPUs per node does not match %d vs. %d"%(num_gpus, self.num_gpus_per_node) - coll_list = self.bte_solver.get_collision_list() - coll_names = self.bte_solver.get_collision_names() if (use_gpu==1): parla_placement = [gpu(gidx_to_device_map(grid_idx,n_grids)) for grid_idx in range(n_grids)] else: parla_placement = [cpu for grid_idx in range(n_grids)] - if csv_write: - data_csv = np.empty((self.tps_npts, 8 + len(coll_list))) + ts = TaskSpace("T") for grid_idx in self.active_grid_idx: @@ -1352,9 +1363,10 @@ def io_output_data(self, grid_idx, u0, plot_data:bool, export_csv:bool, fname:st h_curr = self.bte_solver.normalized_distribution(grid_idx, h_curr) ff = h_curr qoi = self.bte_solver.compute_QoIs(grid_idx, h_curr, effective_mobility=False) - coll_list = self.bte_solver.get_collision_list() - coll_names = self.bte_solver.get_collision_names() - cs_data = self.bte_solver.get_cross_section_data() + + coll_list = self.bte_solver.get_collision_list()[grid_idx] + coll_names = self.bte_solver.get_collision_names()[grid_idx] + cs_data = self.bte_solver.get_cross_section_data()[grid_idx] cs_species = list() for col_idx, (k,v) in enumerate(cs_data.items()): @@ -1781,8 +1793,7 @@ def driver_wo_parla(comm): def __main__(): # TPS solver profile_tt[pp.TPS_SETUP].start() - #tps = libtps.Tps(comm) - tps = libtps.Tps() + tps = libtps.Tps(comm) tps.parseCommandLineArgs(sys.argv) tps.parseInput() tps.chooseDevices() @@ -1829,8 +1840,8 @@ def __main__(): print("tps steps per cycle : ", tps_sper_cycle, "bte_steps per cycle", bte_sper_cycle, flush=True) while (iter Date: Thu, 22 May 2025 09:44:11 -0500 Subject: [PATCH 24/32] Refactor of the python packages --- src/pytps/__init__.py | 11 +++ src/pytps/tabulated_reaction_rates.py | 61 ++++++++++++ src/pytps/utils.py | 5 + src/tps-time-loop.py | 128 ++------------------------ src/tps.py | 12 +-- test/vpath.sh | 5 + 6 files changed, 98 insertions(+), 124 deletions(-) create mode 100644 src/pytps/__init__.py create mode 100644 src/pytps/tabulated_reaction_rates.py create mode 100644 src/pytps/utils.py diff --git a/src/pytps/__init__.py b/src/pytps/__init__.py new file mode 100644 index 000000000..d119d3c3c --- /dev/null +++ b/src/pytps/__init__.py @@ -0,0 +1,11 @@ +from .tabulated_reaction_rates import TabulatedSolver +from .utils import master_print + +# set path to C++ TPS library + +import pathlib +import sys + +path = pathlib.Path(__file__).parent.resolve() +sys.path.append(path + "/.libs") +from libtps import TPS \ No newline at end of file diff --git a/src/pytps/tabulated_reaction_rates.py b/src/pytps/tabulated_reaction_rates.py new file mode 100644 index 000000000..ea360f01f --- /dev/null +++ b/src/pytps/tabulated_reaction_rates.py @@ -0,0 +1,61 @@ +import h5py +import scipy +import scipy.interpolate + +class TabulatedSolver: + def __init__(self, comm, config): + self.comm = comm + self.config = config + self.species_densities = None + self.efield = None + self.heavy_temperature = None + + self.tables = self._read_tables() + self.rates = [] + + def _findPythonReactions(self): + filenames = [] + nreactions = self.config.getint("reactions","number_of_reactions",fallback=0) + for ir in range(nreactions): + sublist = self.config["reactions/reaction{0:d}".format(ir+1)] + rtype = sublist["model"] + if rtype == "bte": + filenames.append(sublist["tabulated/filename"].strip("'")) + + return filenames + + def _read_tables(self): + filenames = self._findPythonReactions() + + tables = [] + for filename in filenames: + with h5py.File(filename, 'r') as fid: + Tcoeff = fid['table'][:] + + tables.append(scipy.interpolate.interp1d(Tcoeff[:,0], Tcoeff[:,1], kind='linear', + bounds_error=False, fill_value='extrapolate')) + + return tables + + def fetch(self, interface): + n_reactions =interface.nComponents(libtps.t2bIndex.ReactionRates) + for r in range(n_reactions): + master_print(self.comm,"Reaction ", r+1, ": ", interface.getReactionEquation(r)) + self.species_densities = np.array(interface.HostRead(libtps.t2bIndex.SpeciesDensities), copy=False) + self.efield = np.array(interface.HostRead(libtps.t2bIndex.ElectricField), copy=False) + self.heavy_temperature = np.array(interface.HostRead(libtps.t2bIndex.HeavyTemperature), copy=False) + + efieldAngularFreq = interface.EfieldAngularFreq() + master_print(self.comm,"Electric field angular frequency: ", efieldAngularFreq) + + def solve(self): + self.rates = [] + for table in self.tables: + self.rates.append(table(self.heavy_temperature)) + + def push(self, interface): + rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False) + offset = 0 + for rate in self.rates: + rates[offset:offset+rate.shape[0]] = rate + offset = offset+rate.shape[0] diff --git a/src/pytps/utils.py b/src/pytps/utils.py new file mode 100644 index 000000000..b46e90563 --- /dev/null +++ b/src/pytps/utils.py @@ -0,0 +1,5 @@ +from mpi4py import MPI + +def master_print(comm: MPI.Comm, *args, **kwargs) -> None: + if comm.rank == 0: + print(*args, **kwargs) diff --git a/src/tps-time-loop.py b/src/tps-time-loop.py index a37c48156..62f87aaba 100755 --- a/src/tps-time-loop.py +++ b/src/tps-time-loop.py @@ -1,129 +1,21 @@ #!/usr/bin/env python3 import sys -import os +import pathlib import numpy as np -import h5py -import scipy -import scipy.interpolate import configparser from mpi4py import MPI -def master_print(comm: MPI.Comm, *args, **kwargs) -> None: - if comm.rank == 0: - print(*args, **kwargs) +# set path to pyTPS library +path = pathlib.Path(__file__).parent.resolve() +sys.path.append(path) +import pytps -class NullSolver: - def __init__(self, comm): - self.comm = comm - self.species_densities = None - self.efield = None - self.heavy_temperature = None - #Reaction 1: 'Ar + E => Ar.+1 + 2 E', - #Reaction 2: 'Ar.+1 + 2 E => Ar + E' - - - def fetch(self, interface): - n_reactions =interface.nComponents(libtps.t2bIndex.ReactionRates) - for r in range(n_reactions): - master_print(self.comm, "Reaction ", r+1, ": ", interface.getReactionEquation(r)) - self.species_densities = np.array(interface.HostRead(libtps.t2bIndex.SpeciesDensities), copy=False) - self.efield = np.array(interface.HostRead(libtps.t2bIndex.ElectricField), copy=False) - self.heavy_temperature = np.array(interface.HostRead(libtps.t2bIndex.HeavyTemperature), copy=False) - - efieldAngularFreq = interface.EfieldAngularFreq() - master_print(self.comm,"Species densities:", self.species_densities.min(), " ", self.species_densities.max()) - master_print(self.comm,"Heavy Temp:", self.heavy_temperature.min(), " ", self.heavy_temperature.max()) - master_print(self.comm,"Efield:", self.efield.min(), " ", self.efield.max()) - master_print(self.comm,"Electric field angular frequency: ", efieldAngularFreq) - - - - def solve(self): - pass - - def push(self, interface): - rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False) - - n_reactions =interface.nComponents(libtps.t2bIndex.ReactionRates) - for r in range(n_reactions): - rates[r*self.heavy_temperature.shape[0]:(r+1)*self.heavy_temperature.shape[0]] = (10.**(-r))*1e-6*self.heavy_temperature - - -class TabulatedSolver: - def __init__(self, comm, config): - self.comm = comm - self.config = config - self.species_densities = None - self.efield = None - self.heavy_temperature = None - - self.tables = self._read_tables() - self.rates = [] - - def _findPythonReactions(self): - filenames = [] - nreactions = self.config.getint("reactions","number_of_reactions",fallback=0) - for ir in range(nreactions): - sublist = self.config["reactions/reaction{0:d}".format(ir+1)] - rtype = sublist["model"] - if rtype == "bte": - filenames.append(sublist["tabulated/filename"].strip("'")) - - return filenames - - def _read_tables(self): - filenames = self._findPythonReactions() - - #["./rate-coefficients/Ionization_Ground.h5", - # "./rate-coefficients/Ionization_Lumped.h5", - # "./rate-coefficients/Excitation_Lumped.h5"] - tables = [] - for filename in filenames: - with h5py.File(filename, 'r') as fid: - Tcoeff = fid['table'][:] - - tables.append(scipy.interpolate.interp1d(Tcoeff[:,0], Tcoeff[:,1], kind='linear', - bounds_error=False, fill_value='extrapolate')) - - return tables - - def fetch(self, interface): - n_reactions =interface.nComponents(libtps.t2bIndex.ReactionRates) - for r in range(n_reactions): - master_print(self.comm,"Reaction ", r+1, ": ", interface.getReactionEquation(r)) - self.species_densities = np.array(interface.HostRead(libtps.t2bIndex.SpeciesDensities), copy=False) - self.efield = np.array(interface.HostRead(libtps.t2bIndex.ElectricField), copy=False) - self.heavy_temperature = np.array(interface.HostRead(libtps.t2bIndex.HeavyTemperature), copy=False) - - efieldAngularFreq = interface.EfieldAngularFreq() - master_print(self.comm,"Electric field angular frequency: ", efieldAngularFreq) - - def solve(self): - self.rates = [] - for table in self.tables: - self.rates.append(table(self.heavy_temperature)) - - def push(self, interface): - rates = np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False) - offset = 0 - for rate in self.rates: - rates[offset:offset+rate.shape[0]] = rate - offset = offset+rate.shape[0] - - - -# 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 = pytps.Tps(comm) tps.parseCommandLineArgs(sys.argv) tps.parseInput() @@ -144,14 +36,14 @@ def push(self, interface): config = configparser.ConfigParser() config.read(ini_name) -boltzmann = TabulatedSolver(comm, config) +boltzmann = pytps.TabulatedSolver(comm, config) -interface = libtps.Tps2Boltzmann(tps) +interface = pytps.Tps2Boltzmann(tps) tps.initInterface(interface) it = 0 max_iters = tps.getRequiredInput("cycle-avg-joule-coupled/max-iters") -master_print(comm,"Max Iters: ", max_iters) +pytps.master_print(comm,"Max Iters: ", max_iters) tps.solveBegin() while it < max_iters: @@ -164,7 +56,7 @@ def push(self, interface): tps.solveStep() it = it+1 - master_print(comm, "it, ", it) + pytps.master_print(comm, "it, ", it) tps.solveEnd() diff --git a/src/tps.py b/src/tps.py index 780bb8c79..1ac38a527 100755 --- a/src/tps.py +++ b/src/tps.py @@ -1,17 +1,17 @@ #!/usr/bin/env python3 +import pathlib import sys -import os from time import perf_counter as time 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 as tps +# set path to pyTPS library +path = pathlib.Path(__file__).parent.resolve() +sys.path.append(path) +import pytps comm = MPI.COMM_WORLD # TPS solver -tps = tps.Tps(comm) +tps = pytps.Tps(comm) tps.parseCommandLineArgs(sys.argv) tps.parseInput() diff --git a/test/vpath.sh b/test/vpath.sh index af1e5eead..7d17260dd 100755 --- a/test/vpath.sh +++ b/test/vpath.sh @@ -29,6 +29,11 @@ if [ ! -d ref_solns ];then ln -s $testDir/ref_solns . fi +# pytps +if [ ! d ../src/pytps ];then + ln -s $testDir/../src/pytps ../src/pytps +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" From 4f3dcf79b0e323fba3313168eb61519f42d68d36 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Thu, 22 May 2025 09:50:11 -0500 Subject: [PATCH 25/32] Fix paths --- src/pytps/__init__.py | 8 +------- src/pytps/tabulated_reaction_rates.py | 3 +++ src/pytps/utils.py | 6 ++++++ 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/pytps/__init__.py b/src/pytps/__init__.py index d119d3c3c..ffce00b21 100644 --- a/src/pytps/__init__.py +++ b/src/pytps/__init__.py @@ -1,11 +1,5 @@ from .tabulated_reaction_rates import TabulatedSolver -from .utils import master_print +from .utils import master_print, libtps -# set path to C++ TPS library -import pathlib -import sys -path = pathlib.Path(__file__).parent.resolve() -sys.path.append(path + "/.libs") -from libtps import TPS \ No newline at end of file diff --git a/src/pytps/tabulated_reaction_rates.py b/src/pytps/tabulated_reaction_rates.py index ea360f01f..04d4f0776 100644 --- a/src/pytps/tabulated_reaction_rates.py +++ b/src/pytps/tabulated_reaction_rates.py @@ -1,6 +1,9 @@ import h5py import scipy import scipy.interpolate +import numpy as np +from .utils import master_print, libtps + class TabulatedSolver: def __init__(self, comm, config): diff --git a/src/pytps/utils.py b/src/pytps/utils.py index b46e90563..b3d73e3b9 100644 --- a/src/pytps/utils.py +++ b/src/pytps/utils.py @@ -1,5 +1,11 @@ from mpi4py import MPI +import pathlib +import sys def master_print(comm: MPI.Comm, *args, **kwargs) -> None: if comm.rank == 0: print(*args, **kwargs) + +path = pathlib.Path(__file__).parent.resolve() +sys.path.append(path + "../.libs") +import libtps From 4f3441c910c575707007d820285708c1deb7e6a6 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Thu, 22 May 2025 09:53:36 -0500 Subject: [PATCH 26/32] Fix structure --- src/tps-time-loop.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tps-time-loop.py b/src/tps-time-loop.py index 62f87aaba..f4f2b4940 100755 --- a/src/tps-time-loop.py +++ b/src/tps-time-loop.py @@ -38,7 +38,7 @@ boltzmann = pytps.TabulatedSolver(comm, config) -interface = pytps.Tps2Boltzmann(tps) +interface = pytps.libtps.Tps2Boltzmann(tps) tps.initInterface(interface) it = 0 From b7281cffef49eabefc887daa05847e515e237534 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Thu, 22 May 2025 12:12:48 -0500 Subject: [PATCH 27/32] Small bugfix --- src/pytps/__init__.py | 2 +- src/pytps/utils.py | 19 ++++++++++++++++--- src/tps-time-loop.py | 17 +++++------------ src/tps.py | 8 +++++--- test/vpath.sh | 2 +- 5 files changed, 28 insertions(+), 20 deletions(-) diff --git a/src/pytps/__init__.py b/src/pytps/__init__.py index ffce00b21..6188ef019 100644 --- a/src/pytps/__init__.py +++ b/src/pytps/__init__.py @@ -1,5 +1,5 @@ from .tabulated_reaction_rates import TabulatedSolver -from .utils import master_print, libtps +from .utils import master_print, resolve_runFile, libtps diff --git a/src/pytps/utils.py b/src/pytps/utils.py index b3d73e3b9..b2e81dbca 100644 --- a/src/pytps/utils.py +++ b/src/pytps/utils.py @@ -1,11 +1,24 @@ from mpi4py import MPI -import pathlib +import os import sys def master_print(comm: MPI.Comm, *args, **kwargs) -> None: if comm.rank == 0: print(*args, **kwargs) -path = pathlib.Path(__file__).parent.resolve() -sys.path.append(path + "../.libs") +def resolve_runFile(cl_args): + ini_name = '' + if '-run' in cl_args: + ini_name = cl_args[cl_args.index('-run') + 1 ] + elif '--runFile' in sys.argv: + ini_name = cl_args[cl_args.index('--runFile') + 1 ] + else: + print("Could not parse command line in python. GOOD BYE!") + exit(-1) + return ini_name + + +base = os.path.abspath(os.path.dirname( sys.argv[0] ) ) +print(os.path.join(base, ".libs")) +sys.path.append(os.path.join(base, ".libs")) import libtps diff --git a/src/tps-time-loop.py b/src/tps-time-loop.py index f4f2b4940..1aa62cd75 100755 --- a/src/tps-time-loop.py +++ b/src/tps-time-loop.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 import sys -import pathlib +import os import numpy as np import configparser @@ -8,14 +8,15 @@ from mpi4py import MPI # set path to pyTPS library -path = pathlib.Path(__file__).parent.resolve() +path = os.path.dirname( os.path.abspath(sys.argv[0]) ) +print(path) sys.path.append(path) import pytps comm = MPI.COMM_WORLD # TPS solver -tps = pytps.Tps(comm) +tps = pytps.libtps.Tps(comm) tps.parseCommandLineArgs(sys.argv) tps.parseInput() @@ -23,15 +24,7 @@ tps.chooseSolver() tps.initialize() -ini_name = '' -if '-run' in sys.argv: - ini_name = sys.argv[sys.argv.index('-run') + 1 ] -elif '--runFile' in sys.argv: - ini_name = sys.argv[sys.argv.index('--runFile') + 1 ] -else: - print("Could not parse command line in python. GOOD BYE!") - exit(-1) - +ini_name = pytps.resolve_runFile(sys.argv) print(ini_name) config = configparser.ConfigParser() config.read(ini_name) diff --git a/src/tps.py b/src/tps.py index 1ac38a527..d6fa7cf69 100755 --- a/src/tps.py +++ b/src/tps.py @@ -1,17 +1,19 @@ #!/usr/bin/env python3 -import pathlib +import os import sys from time import perf_counter as time from mpi4py import MPI # set path to pyTPS library -path = pathlib.Path(__file__).parent.resolve() +path = os.path.dirname( os.path.abspath(sys.argv[0]) ) +print(path) sys.path.append(path) import pytps + comm = MPI.COMM_WORLD # TPS solver -tps = pytps.Tps(comm) +tps = pytps.libtps.Tps(comm) tps.parseCommandLineArgs(sys.argv) tps.parseInput() diff --git a/test/vpath.sh b/test/vpath.sh index 7d17260dd..1f323d5c2 100755 --- a/test/vpath.sh +++ b/test/vpath.sh @@ -30,7 +30,7 @@ if [ ! -d ref_solns ];then fi # pytps -if [ ! d ../src/pytps ];then +if [ ! -d ../src/pytps ];then ln -s $testDir/../src/pytps ../src/pytps fi From c03acc89051013632ca79e1b7500c06cf45ee985 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Thu, 22 May 2025 12:17:28 -0500 Subject: [PATCH 28/32] Docker image --- docker/test/Dockerfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docker/test/Dockerfile b/docker/test/Dockerfile index d22131292..7f3604745 100644 --- a/docker/test/Dockerfile +++ b/docker/test/Dockerfile @@ -142,7 +142,7 @@ RUN . /etc/profile.d/lmod.sh \ RUN rm /tmp/doxygen-1.9.5.src.tar.gz RUN yum -y install graphviz -RUN pip3 install gcovr +RUN . /etc/profile.d/lmod.sh \ && pip3 install gcovr RUN yum -y install gsl-gnu9-ohpc From 077491abddaae53124b248684fd7543bbf540786 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Thu, 22 May 2025 12:35:17 -0500 Subject: [PATCH 29/32] Fix some issues with python imports --- src/tps.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/tps.py b/src/tps.py index d6fa7cf69..68af25b94 100755 --- a/src/tps.py +++ b/src/tps.py @@ -6,14 +6,13 @@ # set path to pyTPS library path = os.path.dirname( os.path.abspath(sys.argv[0]) ) -print(path) -sys.path.append(path) -import pytps +sys.path.append(os.path.join(path, ".libs") +import libtps comm = MPI.COMM_WORLD # TPS solver -tps = pytps.libtps.Tps(comm) +tps = libtps.Tps(comm) tps.parseCommandLineArgs(sys.argv) tps.parseInput() From 760bfb484abf814d591756a5ee08c7e7afc2c4d1 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Thu, 22 May 2025 12:36:59 -0500 Subject: [PATCH 30/32] Fix some issues with python imports --- src/tps.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tps.py b/src/tps.py index 68af25b94..4ca9763e9 100755 --- a/src/tps.py +++ b/src/tps.py @@ -6,7 +6,7 @@ # set path to pyTPS library path = os.path.dirname( os.path.abspath(sys.argv[0]) ) -sys.path.append(os.path.join(path, ".libs") +sys.path.append(os.path.join(path, ".libs")) import libtps From 8c1553c9e9db914de1f9acddf96b385fe05d655a Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Thu, 22 May 2025 13:13:45 -0500 Subject: [PATCH 31/32] Initialize reaction rates from temperature tables --- src/tps-bte_0d3v.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/src/tps-bte_0d3v.py b/src/tps-bte_0d3v.py index 91836a5bf..73f20e6b4 100755 --- a/src/tps-bte_0d3v.py +++ b/src/tps-bte_0d3v.py @@ -58,7 +58,7 @@ def min_mean_max(a, comm: MPI.Comm): # set path to C++ TPS library path = os.path.abspath(os.path.dirname(sys.argv[0])) -sys.path.append(path + "/.libs") +sys.path.append(path) sys.path.append(path + "/../../torch-chemistry") import argon.scripts.synthetic_cs as synthetic_cs sys.path.append(path + "/../../boltzmann/BESolver/python") @@ -66,8 +66,9 @@ def min_mean_max(a, comm: MPI.Comm): crs_folder = path + "/../../torch-chemistry/argon/scripts/crs_lxcat" if not os.path.exists(crs_folder): os.makedirs(crs_folder) - -import libtps + +import pytps +from pytps import libtps, TabulatedSolver from bte_0d3v_batched import bte_0d3v_batched as BoltzmannSolver import utils as bte_utils print(bte_utils) @@ -1804,6 +1805,18 @@ def __main__(): interface = libtps.Tps2Boltzmann(tps) tps.initInterface(interface) tps.solveBegin() + + ## Hack to set-up the correct initial condition + tps.push(interface) + ini_name = pytps.resolve_runFile(sys.argv) + + print(ini_name) + config = configparser.ConfigParser() + config.read(ini_name) + tabulatedSolver = TabulatedSolver(comm, config) + tabulatedSolver.fetch(interface) + tabulatedSolver.solve() + tabulatedSolver.push(interface) # --- first TPS step is needed to initialize the EM fields tps.solveStep() tps.push(interface) From ca3a13eb15b3c57e548e0046cce415c1a32f437a Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Thu, 22 May 2025 13:17:47 -0500 Subject: [PATCH 32/32] Missing fetch --- src/tps-bte_0d3v.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/tps-bte_0d3v.py b/src/tps-bte_0d3v.py index 73f20e6b4..7fd4e8023 100755 --- a/src/tps-bte_0d3v.py +++ b/src/tps-bte_0d3v.py @@ -1817,6 +1817,7 @@ def __main__(): tabulatedSolver.fetch(interface) tabulatedSolver.solve() tabulatedSolver.push(interface) + tps.fetch(interface) # --- first TPS step is needed to initialize the EM fields tps.solveStep() tps.push(interface)