From 7e904573ec91d4107b90074e9ea399bfbccbe71a Mon Sep 17 00:00:00 2001 From: j2jason32 Date: Fri, 24 Oct 2025 14:47:59 -0500 Subject: [PATCH] Surface Response --- doc/source/tutorials/lbs/index.rst | 1 - .../lbs_problem/io/angular_flux_io.cc | 333 +++++++++++++++++- .../lbs_problem/io/lbs_problem_io.h | 49 +++ .../response_evaluator/response_evaluator.cc | 107 +++++- .../response_evaluator/response_evaluator.h | 16 +- python/lib/response.cc | 17 + python/lib/solver.cc | 71 ++++ 7 files changed, 587 insertions(+), 7 deletions(-) diff --git a/doc/source/tutorials/lbs/index.rst b/doc/source/tutorials/lbs/index.rst index d04dd593cb..4b6b5292b2 100644 --- a/doc/source/tutorials/lbs/index.rst +++ b/doc/source/tutorials/lbs/index.rst @@ -1,4 +1,3 @@ - Linear Boltzmann Solver ======================= diff --git a/modules/linear_boltzmann_solvers/lbs_problem/io/angular_flux_io.cc b/modules/linear_boltzmann_solvers/lbs_problem/io/angular_flux_io.cc index edd7872753..b228164773 100644 --- a/modules/linear_boltzmann_solvers/lbs_problem/io/angular_flux_io.cc +++ b/modules/linear_boltzmann_solvers/lbs_problem/io/angular_flux_io.cc @@ -55,7 +55,7 @@ LBSSolverIO::WriteAngularFluxes( const auto& nodes = discretization.GetCellNodeLocations(cell); for (const auto& node : nodes) - { + { nodes_x.push_back(node.x); nodes_y.push_back(node.y); nodes_z.push_back(node.z); @@ -228,9 +228,18 @@ LBSSolverIO::ReadAngularFluxes( H5ReadDataset1D(file_id, group_name + "/values", values); for (uint64_t c = 0; c < file_num_local_cells; ++c) { + // bool isBndry = false; const auto cell_global_id = file_cell_ids[c]; const auto& cell = grid->cells[cell_global_id]; + + const auto& unit_cell_matrices = do_problem.GetUnitCellMatrices(); + const auto& fe_values = unit_cell_matrices.at(cell.local_id); + for (uint64_t i = 0; i < discretization.GetCellNumNodes(cell); ++i) + { + const auto& cell_mapping = discretization.GetCellMapping(cell); + const auto& node_locations = cell_mapping.GetNodeLocations(); + const auto& node_vec = node_locations[i]; for (uint64_t n = 0; n < num_gs_dirs; ++n) for (uint64_t g = 0; g < num_gs_groups; ++g) { @@ -239,9 +248,331 @@ LBSSolverIO::ReadAngularFluxes( psi[dof_map] = values[v]; ++v; } + } + } + } + H5Fclose(file_id); +} + +void +LBSSolverIO::WriteSurfaceAngularFluxes( + DiscreteOrdinatesProblem& do_problem, + const std::string& file_base, + std::vector& bndrys, + std::optional> surfaces) +{ + // Open the HDF5 file + std::string file_name = file_base + std::to_string(opensn::mpi_comm.rank()) + ".h5"; + hid_t file_id = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + OpenSnLogicalErrorIf(file_id < 0, "WriteSurfaceAngularFluxes: Failed to open " + file_name + "."); + + log.Log() << "Writing surface angular flux to " << file_base; + + std::vector>& psi = do_problem.GetPsiNewLocal(); + const auto& supported_bd_ids = do_problem.supported_boundary_ids; + const auto& supported_bd_names = do_problem.supported_boundary_names; + + // Write macro info + const auto& grid = do_problem.GetGrid(); + const auto& discretization = do_problem.GetSpatialDiscretization(); + const auto& groupsets = do_problem.GetGroupsets(); + + auto num_local_cells = grid->local_cells.size(); + auto num_local_nodes = discretization.GetNumLocalNodes(); + auto num_groupsets = groupsets.size(); + + // Store groupsets + H5CreateAttribute(file_id, "num_groupsets", num_groupsets); + + // Check the boundary IDs + std::vector bndry_ids; + const auto unique_bids = grid->GetUniqueBoundaryIDs(); + for (const auto& bndry : bndrys) + { + // Verify if supplied boundary has a valid boundary ID + const auto bndry_id = supported_bd_names.at(bndry); + bndry_ids.push_back(bndry_id); + const auto id = std::find(unique_bids.begin(), unique_bids.end(), bndry_id); + OpenSnInvalidArgumentIf(id == unique_bids.end(), + "Boundary " + bndry + "not found on grid."); + } + + std::vector bndry_tags; + + // =========================================================== + // The goal is limit the number of itrations over all cells + // Sweep through once and map cell_ids to boundaries then + // upack them into 1D vectors for exporting to H5 + // =========================================================== + // The mesh map is structured as follows; + // | face_i | face_i+1 | ... | face_n | + // ----------------------------------------------------------- + // Bndry_id | Cell_i | [nodes_i, nodes_i+1, ..., nodes_n] | + // | Cell_i+1 | [nodes_i, nodes_i+1, ..., nodes_n] | + // | ... | [...] | + // | Cell_n | [...] | + std::map>> mesh_map; + std::map> cell_map, node_map; + std::map> x_map, y_map, z_map; + + // Go through each groupset + unsigned int gset = 0; + for (const auto& groupset : groupsets) + { + // Write groupset info + std::map> surf_map; + + std::map>> data_map; + std::map> mass_map; + + const auto& uk_man = groupset.psi_uk_man_; + const auto& quadrature = groupset.quadrature; + + auto groupset_id = groupset.id; + auto num_gs_dirs = quadrature->omegas.size(); + auto num_gs_groups = groupset.groups.size(); + + // Loop over all cells + for (const auto& cell : grid->local_cells) + { + const uint64_t& cell_id = cell.global_id; + const auto& cell_mapping = discretization.GetCellMapping(cell); + const auto& node_locations = cell_mapping.GetNodeLocations(); + uint64_t num_cell_nodes = 0; + + const auto& unit_cell_matrices = do_problem.GetUnitCellMatrices(); + const auto& fe_values = unit_cell_matrices.at(cell.local_id); + + // Loop over each face of the cell + unsigned int f = 0; + for (const auto& face : cell.faces) + { + bool isSurf = false; + std::string bndry_name; + + // Surface Mapping + if (surfaces.has_value()) + { + const std::string& surf_id = surfaces->first; + double slice = surfaces->second; + + const auto num_face_nodes = cell_mapping.GetNumFaceNodes(f); + for (unsigned int fi = 0; fi < num_face_nodes; ++fi) + { + const auto i = cell_mapping.MapFaceNode(f, fi); + const auto& node_vec = node_locations[i]; + if (node_vec.z == slice) + { + const auto& omega_0 = quadrature->omegas[0]; + const auto mu_0 = omega_0.Dot(face.normal); + bndry_name = surf_id + (mu_0 > 0 ? "_u" : "_d"); + isSurf = true; + bndry_tags.push_back(bndry_name); + } + } + } + // + + // Boundary Mapping + const auto it = std::find(bndry_ids.begin(), bndry_ids.end(), face.neighbor_id); + if (not face.has_neighbor and it != bndry_ids.end()) + { + bndry_name = supported_bd_ids.at(*it); + isSurf = true; + bndry_tags.push_back(bndry_name); + } + + // Write Surface Data + if (isSurf) + { + const auto& int_f_shape_i = fe_values.intS_shapeI[f]; + const auto& M_ij = fe_values.intS_shapeI_shapeJ[f]; + const uint64_t& num_face_nodes = cell_mapping.GetNumFaceNodes(f); + + mesh_map[bndry_name].push_back({cell_id, num_face_nodes}); + cell_map[bndry_name].push_back(cell_id); + node_map[bndry_name].push_back(num_face_nodes); + + num_cell_nodes += num_face_nodes; + for (unsigned int fi = 0; fi < num_face_nodes; ++fi) + { + uint64_t i = cell_mapping.MapFaceNode(f, fi); + const auto& node_vec = node_locations[i]; + if (gset == 0) + { + x_map[bndry_name].push_back(node_vec[0]); + y_map[bndry_name].push_back(node_vec[1]); + z_map[bndry_name].push_back(node_vec[2]); + } + + for (unsigned int d = 0; d < num_gs_dirs; ++d) + { + const auto& omega_d = quadrature->omegas[d]; + const auto weight_d = quadrature->weights[d]; + const auto mu_d = omega_d.Dot(face.normal); + + std::vector data_vec; + data_vec.insert(data_vec.end(), {omega_d.x, omega_d.y, omega_d.z}); + data_vec.push_back(mu_d); + data_vec.push_back(weight_d); + data_vec.push_back(int_f_shape_i(i)); + for (uint64_t g = 0; g < num_gs_groups; ++g) + { + const auto dof_map = discretization.MapDOFLocal(cell, i, uk_man, d, g); + data_vec.push_back(psi[groupset_id][dof_map]); + } + // Move the vector to avoid unecessary copy + data_map[bndry_name].push_back(std::move(data_vec)); + } + + for (unsigned int fj = 0; fj < num_face_nodes; ++fj) + { + const auto j = cell_mapping.MapFaceNode(f, fj); + mass_map[bndry_name].push_back(M_ij(i, j)); + } + } + } + ++f; + } + } + + // Export data to HDF5 + std::string group_name = "groupset_" + std::to_string(groupset_id); + H5CreateGroup(file_id, group_name); + H5CreateAttribute(file_id, group_name + "/num_directions", num_gs_dirs); + H5CreateAttribute(file_id, group_name + "/num_groups", num_gs_groups); + + H5CreateGroup(file_id, "mesh"); + for (const auto& bndry_id : bndry_tags) + { + if (gset == 0) + { + std::cout << "Boundry ID:" << bndry_id << std::endl; + const auto& cell_ids = cell_map.at(bndry_id); + const auto& num_nodes = node_map.at(bndry_id); + const auto& x_bndry = x_map.at(bndry_id); + const auto& y_bndry = y_map.at(bndry_id); + const auto& z_bndry = z_map.at(bndry_id); + + std::string bndry_mesh = std::string("mesh/") + bndry_id; + H5CreateGroup(file_id, bndry_mesh); + H5WriteDataset1D(file_id, bndry_mesh + "/cell_ids", cell_ids); + H5WriteDataset1D(file_id, bndry_mesh + "/num_nodes", num_nodes); + H5WriteDataset1D(file_id, bndry_mesh + "/nodes_x", x_bndry); + H5WriteDataset1D(file_id, bndry_mesh + "/nodes_y", y_bndry); + H5WriteDataset1D(file_id, bndry_mesh + "/nodes_z", z_bndry); + } + + std::string bndry_grp = group_name + std::string("/") + bndry_id; + H5CreateGroup(file_id, bndry_grp); + + // Write the groupset surface angular flux data + std::vector omega; + std::vector mu; + std::vector wt_d; + std::vector fe_shape; + std::vector surf_flux; + + const auto& data_vectors = data_map.at(bndry_id); + for (const auto&vec : data_vectors) + { + omega.insert(omega.end(), {vec[0], vec[1], vec[2]}); + mu.push_back(vec[3]); + wt_d.push_back(vec[4]); + fe_shape.push_back(vec[5]); + surf_flux.insert(surf_flux.end(), vec.begin()+6, vec.end()); + } + H5WriteDataset1D(file_id, bndry_grp + "/omega", omega); + H5WriteDataset1D(file_id, bndry_grp + "/mu", mu); + H5WriteDataset1D(file_id, bndry_grp + "/wt_d", wt_d); + H5WriteDataset1D(file_id, bndry_grp + "/fe_shape", fe_shape); + H5WriteDataset1D(file_id, bndry_grp + "/surf_flux", surf_flux); + + // Write mass matrix information + const auto& mass_vector = mass_map.at(bndry_id); + H5WriteDataset1D(file_id, bndry_grp + "/M_ij", mass_vector); } + ++gset; } + + ssize_t num_open_objs = H5Fget_obj_count(file_id, H5F_OBJ_ALL); H5Fclose(file_id); } +std::vector +LBSSolverIO::ReadSurfaceAngularFluxes( + DiscreteOrdinatesProblem& do_problem, + const std::string& file_base, + std::vector& bndrys) +{ + std::vector surf_fluxes; + + // Open HDF5 file + std::string file_name = file_base + std::to_string(opensn::mpi_comm.rank()) + ".h5"; + hid_t file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + OpenSnLogicalErrorIf(file_id < 0, "Failed to open " + file_name + "."); + + log.Log() << "Reading surface angular flux file from " << file_base; + + const auto& supported_bd_ids = do_problem.supported_boundary_ids; + const auto& supported_bd_names = do_problem.supported_boundary_names; + + // Read macro data and check for compatibility + uint64_t file_num_groupsets; + H5ReadAttribute(file_id, "num_groupsets", file_num_groupsets); + + const auto& grid = do_problem.GetGrid(); + const auto& discretization = do_problem.GetSpatialDiscretization(); + const auto& groupsets = do_problem.GetGroupsets(); + auto num_groupsets = groupsets.size(); + + OpenSnLogicalErrorIf(file_num_groupsets != num_groupsets, + "Incompatible number of groupsets found in file " + file_name + "."); + + // Go through each groupset + for (const auto& groupset : groupsets) + { + const auto& uk_man = groupset.psi_uk_man_; + const auto& quadrature = groupset.quadrature; + + auto groupset_id = groupset.id; + auto num_gs_dirs = quadrature->omegas.size(); + auto num_gs_groups = groupset.groups.size(); + + uint64_t file_num_gs_dirs; + uint64_t file_num_gs_groups; + std::string group_name = "groupset_" + std::to_string(groupset_id); + H5ReadAttribute(file_id, group_name + "/num_directions", file_num_gs_dirs); + H5ReadAttribute(file_id, group_name + "/num_groups", file_num_gs_groups); + + for (const auto& bndry : bndrys) + { + SurfaceMap surf_map; + SurfaceData surf_data; + SurfaceAngularFlux surf_flux; + + std::cout << "Reading Surface: " << bndry << std::endl; + std::string mesh_tag = "mesh/" + bndry; + H5ReadDataset1D(file_id, mesh_tag + "/cell_ids", surf_map.cell_ids); + H5ReadDataset1D(file_id, mesh_tag + "/num_nodes", surf_map.num_nodes); + H5ReadDataset1D(file_id, mesh_tag + "/nodes_x", surf_map.nodes_x); + H5ReadDataset1D(file_id, mesh_tag + "/nodes_y", surf_map.nodes_y); + H5ReadDataset1D(file_id, mesh_tag + "/nodes_z", surf_map.nodes_z); + std::string bndry_grp = group_name + "/" + bndry; + H5ReadDataset1D(file_id, bndry_grp + "/omega", surf_data.omega); + H5ReadDataset1D(file_id, bndry_grp + "/mu", surf_data.mu); + H5ReadDataset1D(file_id, bndry_grp + "/wt_d", surf_data.wt_d); + H5ReadDataset1D(file_id, bndry_grp + "/M_ij", surf_data.M_ij); + H5ReadDataset1D(file_id, bndry_grp + "/surf_flux", surf_data.psi); + + surf_flux.mapping = std::move(surf_map); + surf_flux.data = std::move(surf_data); + surf_fluxes.push_back(std::move(surf_flux)); + } + } + H5Fclose(file_id); + + return surf_fluxes; +} + } // namespace opensn diff --git a/modules/linear_boltzmann_solvers/lbs_problem/io/lbs_problem_io.h b/modules/linear_boltzmann_solvers/lbs_problem/io/lbs_problem_io.h index d55b0f7043..cae5a3bbaa 100644 --- a/modules/linear_boltzmann_solvers/lbs_problem/io/lbs_problem_io.h +++ b/modules/linear_boltzmann_solvers/lbs_problem/io/lbs_problem_io.h @@ -3,9 +3,11 @@ #pragma once +#include #include #include #include +#include #include namespace opensn @@ -43,6 +45,53 @@ class LBSSolverIO std::optional>>> opt_dest = std::nullopt); + struct SurfaceMap { + std::vector cell_ids; + std::vector num_nodes; + std::vector nodes_x; + std::vector nodes_y; + std::vector nodes_z; + }; + + struct SurfaceData { + std::vector omega; + std::vector mu; + std::vector wt_d; + std::vector M_ij; + std::vector psi; + }; + + struct SurfaceAngularFlux { + SurfaceMap mapping; + SurfaceData data; + }; + + /** + * Write surface angular flux vector(s) to a file. + * + * \param lbs_solver LBS solver + * \param file_base File name stem + * \param bndrys Map of boundary names and ids + */ + static void WriteSurfaceAngularFluxes( + DiscreteOrdinatesProblem& do_problem, + const std::string& file_stem, + std::vector& bndrys, + std::optional> surfaces); + + /** + * Read a surface angular flux vector from a file. + * + * \param lbs_problem LBS problem + * \param file_base File name stem + * \param per_material Optional angular flux destination vector + */ + static std::vector + ReadSurfaceAngularFluxes( + DiscreteOrdinatesProblem& do_problem, + const std::string& file_stem, + std::vector& bndrys); + /** * Write a flux moments vector to a file. * diff --git a/modules/linear_boltzmann_solvers/response_evaluator/response_evaluator.cc b/modules/linear_boltzmann_solvers/response_evaluator/response_evaluator.cc index ec8f08a7d1..a25c56f81f 100644 --- a/modules/linear_boltzmann_solvers/response_evaluator/response_evaluator.cc +++ b/modules/linear_boltzmann_solvers/response_evaluator/response_evaluator.cc @@ -118,6 +118,16 @@ ResponseEvaluator::GetBufferOptionsBlock() "file_prefixes", "A table containing file prefixes for flux moments and angular flux binary files. " "These are keyed by \"flux_moments\" and \"angular_fluxes\", respectively."); + params.AddOptionalParameter( + "bndry", + {}, + "Optional list of boundary ids." + ); + params.AddOptionalParameter( + "surf", + {}, + "Optional list of surface ids." + ); return params; } @@ -144,7 +154,29 @@ ResponseEvaluator::SetBufferOptions(const InputParameters& input) LBSSolverIO::ReadAngularFluxes( *do_problem_, prefixes.GetParamValue("angular_fluxes"), psi); - adjoint_buffers_[name] = {phi, psi}; + std::vector surf_psi; + if (prefixes.Has("surface_angular_fluxes")) + { + const auto bndry = params.GetParamValue("bndry"); + const auto surf = params.GetParamValue("surf"); + + std::vector bndrys; + if (!bndry.empty()) + bndrys.push_back(surf); + else if (!surf.empty()) + { + const std::string surf_u = surf + "_u"; + const std::string surf_d = surf + "_d"; + bndrys.insert(bndrys.end(), {surf_u, surf_d}); + } + + // bndrys.push_back(surf); + + surf_psi = LBSSolverIO::ReadSurfaceAngularFluxes( + *do_problem_, prefixes.GetParamValue("surface_angular_fluxes"), bndrys); + } + + adjoint_buffers_[name] = {phi, psi, surf_psi}; log.Log0Verbose1() << "Adjoint buffer " << name << " added to the stack."; } @@ -315,8 +347,11 @@ double ResponseEvaluator::EvaluateResponse(const std::string& buffer) const { const auto& buffer_data = adjoint_buffers_.at(buffer); - const auto& phi_dagger = buffer_data.first; - const auto& psi_dagger = buffer_data.second; + // const auto& phi_dagger = buffer_data.first; + // const auto& psi_dagger = buffer_data.second; + + const auto& phi_dagger = buffer_data.flux_moments; + const auto& psi_dagger = buffer_data.angular_fluxes; OpenSnLogicalErrorIf(not material_sources_.empty() and phi_dagger.empty(), "If material sources are present, adjoint flux moments " @@ -463,6 +498,72 @@ ResponseEvaluator::EvaluateResponse(const std::string& buffer) const return global_response; } +double +ResponseEvaluator::EvaluateSurfaceResponse(const std::string& fwd_buffer, + const std::string& adj_buffer) const +{ + const auto& fwd_data = adjoint_buffers_.at(fwd_buffer); + const auto& fwd_surfaces = fwd_data.surface_angular_fluxes; + OpenSnLogicalErrorIf(fwd_surfaces.empty(), + "Surface flux data must be available " + "for a surface response evaluation."); + + const auto& adj_data = adjoint_buffers_.at(adj_buffer); + const auto& adj_surfaces = adj_data.surface_angular_fluxes; + OpenSnLogicalErrorIf(adj_surfaces.empty(), + "Surface adjoint flux data must be available " + "for a surface response evaluation."); + + double response = 0.0; + const auto& num_surfs = fwd_surfaces.size() - 1; + for (size_t si=0; si < num_surfs; ++si) + { + const auto& adj = adj_surfaces[1]; + const auto& fwd = fwd_surfaces[0]; + + size_t gi = 0; // start of global index + for (size_t ci=0; ci < fwd.mapping.cell_ids.size(); ++ci) + { + size_t mi = 0; // start of mass index + const auto& num_nodes = fwd.mapping.num_nodes[ci]; + for (size_t ni=0; ni < num_nodes; ++ni) + { + const auto& node_strd = gi + ni; + for (const auto& groupset : do_problem_->GetGroupsets()) + { + auto groupset_id = groupset.id; + auto num_gs_groups = groupset.groups.size(); + + const auto& quadrature = groupset.quadrature; + auto num_gs_dirs = quadrature->omegas.size(); + for (unsigned int d = 0; d < num_gs_dirs; ++d) + { + const auto& dir_strd = node_strd * num_gs_dirs + d; + const auto& adir_strd = node_strd * num_gs_dirs + num_gs_dirs - d - 1; + + const auto& mu_d = adj.data.mu[adir_strd]; + const auto& wt_d = adj.data.wt_d[adir_strd]; + for (uint64_t g = 0; g < num_gs_groups; ++g) + { + const auto& grp_strd = dir_strd * num_gs_groups + g; + const auto& agrp_strd = adir_strd * num_gs_groups + g; + for (size_t nj=0; nj < num_nodes; ++nj) + { + const auto& mass_strd = mi + nj; + response += mu_d * wt_d * adj.data.M_ij[mass_strd] + * adj.data.psi[agrp_strd] * fwd.data.psi[grp_strd]; + } + } + } + } + mi += num_nodes; + } + gi += num_nodes; + } + } + return response; +} + std::vector ResponseEvaluator::EvaluateBoundaryCondition(const uint64_t boundary_id, const Vector3& node, diff --git a/modules/linear_boltzmann_solvers/response_evaluator/response_evaluator.h b/modules/linear_boltzmann_solvers/response_evaluator/response_evaluator.h index 7c7c215f71..022aefab05 100644 --- a/modules/linear_boltzmann_solvers/response_evaluator/response_evaluator.h +++ b/modules/linear_boltzmann_solvers/response_evaluator/response_evaluator.h @@ -5,6 +5,7 @@ #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h" #include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/io/lbs_problem_io.h" #include "framework/parameters/input_parameters.h" #include @@ -58,8 +59,13 @@ class ResponseEvaluator private: using FluxMomentBuffer = std::vector; using AngularFluxBuffer = std::vector>; - using AdjointBuffer = std::pair; - + using SurfaceAngularFluxBuffer = std::vector; + struct AdjointBuffer + { + FluxMomentBuffer flux_moments; + AngularFluxBuffer angular_fluxes; + SurfaceAngularFluxBuffer surface_angular_fluxes; + }; using MaterialSources = std::map>; using PointSources = std::vector>; using VolumetricSources = std::vector>; @@ -95,6 +101,12 @@ class ResponseEvaluator */ double EvaluateResponse(const std::string& buffer_name) const; + /** + * Evaluate a response using the specified forward buffer and + * adjoint buffer towards a specified surface id. + */ + double EvaluateSurfaceResponse(const std::string& fwd_buffer, + const std::string& adj_buffer) const; private: /** * Evaluates a boundary source and returns the angular flux on the boundary. diff --git a/python/lib/response.cc b/python/lib/response.cc index 852f7aa18c..082df5732c 100644 --- a/python/lib/response.cc +++ b/python/lib/response.cc @@ -61,6 +61,23 @@ WrapResEval(py::module& response) )", py::arg("buffer_name") ); + res_eval.def( + "EvaluateSurfaceResponse", + &ResponseEvaluator::EvaluateSurfaceResponse, + R"( + Evaluate a surface response using the specified bndry name and adjoint buffer with + the currently defined forward sources in the solver. + + Parameters + ---------- + fwd_buffer: str + Foward buffer name + adj_buffer: str + Adjoint buffer name + )", + py::arg("fwd_buffer"), + py::arg("adj_buffer") + ); res_eval.def( "SetOptions", [](ResponseEvaluator& self, py::kwargs& params) diff --git a/python/lib/solver.cc b/python/lib/solver.cc index 581766c638..772cc324dd 100644 --- a/python/lib/solver.cc +++ b/python/lib/solver.cc @@ -403,6 +403,77 @@ WrapLBS(py::module& slv) )", py::arg("file_base") ); + lbs_problem.def( + "WriteSurfaceAngularFluxes", + [](DiscreteOrdinatesProblem& self, + const std::string& file_base, + py::list bndry_names, + py::object surfaces) + { + // Map boundary names + std::map allowed_bd_names = LBSProblem::supported_boundary_names; + std::map allowed_bd_ids = LBSProblem::supported_boundary_ids; + std::vector bndrys; + for (py::handle name : bndry_names) + bndrys.push_back(name.cast()); + + // Map surface names + std::optional> opt_surfaces; + if (!surfaces.is_none()) + { + py::list surf_list = surfaces.cast(); + if (py::len(surf_list) == 2) + { + std::string surf_id = surf_list[0].cast(); + double slice = surf_list[1].cast(); + opt_surfaces = std::make_pair(surf_id, slice); + } + } + + // LBSSolverIO::WriteSurfaceAngularFluxes(self, file_base, bndry_map); + LBSSolverIO::WriteSurfaceAngularFluxes(self, file_base, bndrys, opt_surfaces); + }, + R"( + Write surface angular flux data to file. + + Parameters + ---------- + file_base: str + File basename. + )", + py::arg("file_base"), + py::arg("bndry_names"), + py::arg("surfaces") = py::none() + ); + lbs_problem.def( + "ReadSurfaceAngularFluxes", + [](DiscreteOrdinatesProblem& self, const std::string& file_base, py::list bndry_names) + { + std::map supported_bd_names = LBSProblem::supported_boundary_names; + std::map supported_bd_ids = LBSProblem::supported_boundary_ids; + + // std::map bndry_map; + std::vector bndrys; + for (py::handle name : bndry_names) + { + std::string bndry = name.cast(); + // bndry_map[bndry] = supported_bd_names.at(bndry); + bndrys.push_back(bndry); + } + // LBSSolverIO::ReadSurfaceAngularFluxes(self, file_base, bndry_map); + LBSSolverIO::ReadSurfaceAngularFluxes(self, file_base, bndrys); + }, + R"( + Read surface angular fluxes from file. + + Parameters + ---------- + file_base: str + File basename. + )", + py::arg("file_base"), + py::arg("bndry_names") + ); lbs_problem.def( "SetPointSources", [](LBSProblem& self, py::kwargs& params)