25654184358

Committed 09 May 2026 03:58AM UTC coverage: 75.687%. Remained the same

Build # 25654184358

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #1042 from wdhawkins/ang_quadrature_refactor

Refactoring angular quadrature classes.

Coverage Stats

230 of 278 new or added lines in 27 files covered. (82.73%)

181 existing lines in 12 files now uncovered.

22258 of 29408 relevant lines covered (75.69%)

64739266.89 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

54.93

/modules/linear_boltzmann_solvers/response_evaluator/response_evaluator.cc

// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
// SPDX-License-Identifier: MIT

#include "modules/linear_boltzmann_solvers/response_evaluator/response_evaluator.h"
#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h"
#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/io/discrete_ordinates_problem_io.h"
#include "modules/linear_boltzmann_solvers/lbs_problem/point_source/point_source.h"
#include "modules/linear_boltzmann_solvers/lbs_problem/volumetric_source/volumetric_source.h"
#include "modules/linear_boltzmann_solvers/lbs_problem/io/lbs_problem_io.h"
#include "framework/mesh/mesh_continuum/mesh_continuum.h"
#include "framework/logging/log.h"
#include "framework/object_factory.h"
#include "framework/runtime.h"
#include "mpicpp-lite/mpicpp-lite.h"

namespace mpi = mpicpp_lite;

namespace opensn
{

OpenSnRegisterObjectInNamespace(lbs, ResponseEvaluator);

InputParameters
ResponseEvaluator::GetInputParameters()
{
  InputParameters params;
  params.SetGeneralDescription(
    "A utility class for evaluating responses using precomputed adjoint solutions "
    "and arbitrary forward sources.");

  params.AddRequiredParameter<std::shared_ptr<Problem>>("problem",
                                                        "A handle to an existing LBS problem.");
  params.AddOptionalParameterBlock(
    "options", ParameterBlock(), "The specification of adjoint buffers and forward to use.");
  params.LinkParameterToBlock("options", "response::OptionsBlock");

  return params;
}

std::shared_ptr<ResponseEvaluator>
ResponseEvaluator::Create(const ParameterBlock& params)
{
  auto& factory = opensn::ObjectFactory::GetInstance();
  return factory.Create<ResponseEvaluator>("lbs::ResponseEvaluator", params);
}

ResponseEvaluator::ResponseEvaluator(const InputParameters& params)
  : do_problem_(params.GetSharedPtrParam<Problem, DiscreteOrdinatesProblem>("problem"))
{
  if (params.IsParameterValid("options"))
  {
    auto options = GetOptionsBlock();
    options.AssignParameters(params.GetParam("options"));
    SetOptions(options);
  }
}

InputParameters
ResponseEvaluator::GetOptionsBlock()
{
  InputParameters params;
  params.SetGeneralDescription("A block of options for the response evaluator for adding adjoint "
                               "buffers and defining forward sources.");

  params.AddOptionalParameterArray(
    "buffers", {}, "An array of tables containing adjoint buffer specifications.");
  params.LinkParameterToBlock("buffers", "response::BufferOptionsBlock");

  params.AddOptionalParameter("clear_sources", false, "A flag to clear existing sources.");
  params.AddOptionalParameterBlock(
    "sources", ParameterBlock(), "An array of tables containing source specification information.");
  params.LinkParameterToBlock("sources", "response::SourceOptionsBlock");

  return params;
}

void
ResponseEvaluator::SetOptions(const InputParameters& params)
{
  if (params.IsParameterValid("buffers"))
  {
    const auto& user_buffer_params = params.GetParam("buffers");
    user_buffer_params.RequireBlockTypeIs(ParameterBlockType::ARRAY);
    for (int p = 0; p < user_buffer_params.GetNumParameters(); ++p)
    {
      auto buffer_params = GetBufferOptionsBlock();
      buffer_params.AssignParameters(user_buffer_params.GetParam(p));
      SetBufferOptions(buffer_params);
    }
  }

  if (params.IsParameterValid("clear_sources") and params.GetParamValue<bool>("clear_sources"))
  {
    material_sources_.clear();
    point_sources_.clear();
    volumetric_sources_.clear();
    boundary_sources_.clear();
  }

  if (params.IsParameterValid("sources"))
  {
    auto source_params = GetSourceOptionsBlock();
    source_params.AssignParameters(params.GetParam("sources"));
    SetSourceOptions(source_params);
  }
}

InputParameters
ResponseEvaluator::GetBufferOptionsBlock()
{
  InputParameters params;
  params.SetGeneralDescription("Options for adding adjoint buffers to the response evaluator.");

  params.AddRequiredParameter<std::string>(
    "name",
    "A name given to the buffer to identify it when querying the response evaluation routine.");
  params.AddRequiredParameterBlock(
    "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.");

  return params;
}

void
ResponseEvaluator::SetBufferOptions(const InputParameters& input)
{
  auto params = opensn::ResponseEvaluator::GetBufferOptionsBlock();
  params.AssignParameters(input);

  const auto name = params.GetParamValue<std::string>("name");
  OpenSnInvalidArgumentIf(adjoint_buffers_.count(name) > 0,
                          "An adjoint buffer with name " + name + " already exists.");

  const auto prefixes = params.GetParam("file_prefixes");

  std::vector<double> phi;
  if (prefixes.Has("flux_moments"))
    LBSSolverIO::ReadFluxMoments(
      *do_problem_, prefixes.GetParamValue<std::string>("flux_moments"), false, phi);

  std::vector<std::vector<double>> psi;
  if (prefixes.Has("angular_fluxes"))
    DiscreteOrdinatesProblemIO::ReadAngularFluxes(
      *do_problem_, prefixes.GetParamValue<std::string>("angular_fluxes"), psi);

  adjoint_buffers_[name] = {phi, psi};
  log.Log0Verbose1() << "Adjoint buffer " << name << " added to the stack.";
}

InputParameters
ResponseEvaluator::GetSourceOptionsBlock()
{
  InputParameters params;
  params.SetGeneralDescription("A table of various forward source specifications.");

  params.AddOptionalParameterArray(
    "material", {}, "An array of tables containing material source specifications.");
  params.LinkParameterToBlock("material", "response::MaterialSourceOptionsBlock");

  params.AddOptionalParameterArray<std::shared_ptr<PointSource>>(
    "point", {}, "An array of tables containing point source handles.");

  params.AddOptionalParameterArray<std::shared_ptr<VolumetricSource>>(
    "volumetric", {}, "An array of tables containing volumetric source handles.");

  params.AddOptionalParameterArray(
    "boundary", {}, "An array of tables containing boundary source specifications.");
  params.LinkParameterToBlock("boundary", "response::BoundarySourceOptionsBlock");

  return params;
}

InputParameters
ResponseEvaluator::GetBoundarySourceOptionsBlock()
{
  InputParameters params;

  params.SetGeneralDescription("Boundary source specification for response evaluations.");
  params.SetClassName("Response Boundary Source");
  params.AddRequiredParameter<std::string>("name",
                                           "Boundary name that identifies the specific boundary");
  params.AddRequiredParameter<std::string>("type", "Boundary type specification.");
  params.AddOptionalParameterArray<double>("group_strength",
                                           {},
                                           "Required only if \"type\" is \"isotropic\". An array "
                                           "of isotropic strength per group");
  params.ConstrainParameterRange("type", AllowableRangeList::New({"isotropic"}));

  return params;
}

void
ResponseEvaluator::SetSourceOptions(const InputParameters& input)
{
  auto params = ResponseEvaluator::GetSourceOptionsBlock();
  params.AssignParameters(input);

  params.RequireBlockTypeIs(ParameterBlockType::BLOCK);

  // Add material sources
  if (params.Has("material"))
  {
    const auto& user_msrc_params = params.GetParam("material");
    for (int p = 0; p < user_msrc_params.GetNumParameters(); ++p)
    {
      auto msrc_params = GetMaterialSourceOptionsBlock();
      msrc_params.AssignParameters(user_msrc_params.GetParam(p));
      SetMaterialSourceOptions(msrc_params);
    }
  }

  // Add point sources
  if (params.Has("point"))
  {
    const auto& user_psrc_params = params.GetParam("point");
    for (int p = 0; p < user_psrc_params.GetNumParameters(); ++p)
    {
      point_sources_.push_back(
        user_psrc_params.GetParam(p).GetValue<std::shared_ptr<PointSource>>());
      point_sources_.back()->Initialize(*do_problem_);
    }
  }

  // Add volumetric sources
  if (params.Has("volumetric"))
  {
    const auto& user_dsrc_params = params.GetParam("volumetric");
    for (int p = 0; p < user_dsrc_params.GetNumParameters(); ++p)
    {
      volumetric_sources_.push_back(
        user_dsrc_params.GetParam(p).GetValue<std::shared_ptr<VolumetricSource>>());
      volumetric_sources_.back()->Initialize(*do_problem_);
    }
  }

  // Add boundary sources
  if (params.Has("boundary"))
  {
    const auto& user_bsrc_params = params.GetParam("boundary");
    for (int p = 0; p < user_bsrc_params.GetNumParameters(); ++p)
    {
      auto bsrc_params = ResponseEvaluator::GetBoundarySourceOptionsBlock();
      bsrc_params.AssignParameters(user_bsrc_params.GetParam(p));
      SetBoundarySourceOptions(bsrc_params);
    }
  }
}

InputParameters
ResponseEvaluator::GetMaterialSourceOptionsBlock()
{
  InputParameters params;
  params.SetGeneralDescription(
    "Options for adding material-based forward sources to the response evaluator.");

  params.AddRequiredParameter<unsigned int>("block_id", "The block id the source belongs to.");
  params.AddRequiredParameterArray("strength", "The group-wise material source strength.");

  return params;
}

void
ResponseEvaluator::SetMaterialSourceOptions(const InputParameters& params)
{
  const auto blkid = params.GetParamValue<unsigned int>("block_id");
  OpenSnInvalidArgumentIf(material_sources_.count(blkid) > 0,
                          "A material source for block id " + std::to_string(blkid) +
                            " already exists.");

  const auto values = params.GetParamVectorValue<double>("strength");
  OpenSnInvalidArgumentIf(values.size() != do_problem_->GetNumGroups(),
                          "The number of material source values and groups "
                          "in the underlying solver do not match. "
                          "Expected " +
                            std::to_string(do_problem_->GetNumGroups()) + " but got " +
                            std::to_string(values.size()) + ".");

  material_sources_[blkid] = values;
  log.Log0Verbose1() << "Material source for block id " << blkid << " added to the stack.";
}

void
ResponseEvaluator::SetBoundarySourceOptions(const InputParameters& params)
{
  const auto bndry_name = params.GetParamValue<std::string>("name");
  if (params.IsParameterValid("function"))
    throw std::runtime_error("Boundary '" + bndry_name +
                             "' in ResponseEvaluator does not support \"function\".");
  const auto bndry_type = params.GetParamValue<std::string>("type");

  auto grid = do_problem_->GetGrid();
  const auto bnd_name_map = grid->GetBoundaryNameMap();
  OpenSnInvalidArgumentIf(not bnd_name_map.count(bndry_name),
                          "Boundary \"" + bndry_name + "\" does not exist.");
  const auto bid = bnd_name_map.at(bndry_name);
  if (bndry_type == "isotropic")
  {
    OpenSnInvalidArgumentIf(not params.Has("group_strength"),
                            "Parameter \"group_strength\" is required for "
                            "boundaries of type \"isotropic\".");
    params.RequireParameterBlockTypeIs("group_strength", ParameterBlockType::ARRAY);
    const auto values = params.GetParamVectorValue<double>("group_strength");
    OpenSnInvalidArgumentIf(values.size() != do_problem_->GetNumGroups(),
                            "The number of boundary source values does not match the "
                            "number of groups. Expected " +
                              std::to_string(do_problem_->GetNumGroups()) + " but got " +
                              std::to_string(values.size()) + ".");

    boundary_sources_[bid] = {LBSBoundaryType::ISOTROPIC, values};
  }
  else
    log.Log0Warning() << "Unsupported boundary type. Skipping the entry.";
}

void
ResponseEvaluator::ClearForwardSources()
{
  material_sources_.clear();
  point_sources_.clear();
  volumetric_sources_.clear();
  boundary_sources_.clear();
}

void
ResponseEvaluator::AddResponseBuffers(const InputParameters& params)
{
  params.RequireBlockTypeIs(ParameterBlockType::ARRAY);
  for (size_t p = 0; p < params.GetNumParameters(); ++p)
  {
    auto spec = ResponseEvaluator::GetBufferOptionsBlock();
    spec.AssignParameters(params.GetParam(p));
    SetBufferOptions(spec);
  }
}

void
ResponseEvaluator::AddResponseSources(const InputParameters& params)
{
  auto spec = ResponseEvaluator::GetSourceOptionsBlock();
  spec.AssignParameters(params);
  SetSourceOptions(spec);
}

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;

  OpenSnLogicalErrorIf(not material_sources_.empty() and phi_dagger.empty(),
                       "If material sources are present, adjoint flux moments "
                       "must be available for response evaluation.");
  OpenSnLogicalErrorIf(not point_sources_.empty() and phi_dagger.empty(),
                       "If point sources are set, adjoint flux moments "
                       "must be available for response evaluation.");
  OpenSnLogicalErrorIf(not volumetric_sources_.empty() and phi_dagger.empty(),
                       "if volumetric sources are set, adjoint flux moments "
                       "must be available for response evaluation.");
  OpenSnLogicalErrorIf(not boundary_sources_.empty() and psi_dagger.empty(),
                       "If boundary sources are set, adjoint angular fluxes "
                       "must be available for response evaluation.");

  const auto& grid = do_problem_->GetGrid();
  const auto& discretization = do_problem_->GetSpatialDiscretization();
  const auto& transport_views = do_problem_->GetCellTransportViews();
  const auto& unit_cell_matrices = do_problem_->GetUnitCellMatrices();
  const auto num_groups = do_problem_->GetNumGroups();

  double local_response = 0.0;

  // Material sources
  if (not material_sources_.empty())
  {
    for (const auto& cell : grid->local_cells)
    {
      const auto& cell_mapping = discretization.GetCellMapping(cell);
      const auto& transport_view = transport_views[cell.local_id];
      const auto& fe_values = unit_cell_matrices[cell.local_id];
      const auto num_cell_nodes = cell_mapping.GetNumNodes();

      if (material_sources_.count(cell.block_id) > 0)
      {
        const auto& src = material_sources_.at(cell.block_id);
        for (size_t i = 0; i < num_cell_nodes; ++i)
        {
          const auto dof_map = transport_view.MapDOF(i, 0, 0);
          const auto& V_i = fe_values.intV_shapeI(i);
          for (unsigned int g = 0; g < num_groups; ++g)
            local_response += src[g] * phi_dagger[dof_map + g] * V_i;
        }
      }
    } // for cell
  } // if material sources

  // Boundary sources
  if (not boundary_sources_.empty())
  {
    size_t gs = 0;
    for (const auto& groupset : do_problem_->GetGroupsets())
    {
      const auto& uk_man = groupset.psi_uk_man_;
      const auto& quadrature = groupset.quadrature;
      const auto num_gs_angles = quadrature->GetNumAngles();
      const auto& num_gs_groups = groupset.GetNumGroups();

      for (const auto& cell : grid->local_cells)
      {
        const auto& cell_mapping = discretization.GetCellMapping(cell);
        const auto& fe_values = unit_cell_matrices[cell.local_id];

        size_t f = 0;
        for (const auto& face : cell.faces)
        {
          if (not face.has_neighbor and boundary_sources_.count(face.neighbor_id) > 0)
          {
            const auto bndry_id = face.neighbor_id;
            const auto num_face_nodes = cell_mapping.GetNumFaceNodes(f);
            for (size_t fi = 0; fi < num_face_nodes; ++fi)
            {
              const auto i = cell_mapping.MapFaceNode(f, fi);
              const auto& node = grid->vertices[cell.vertex_ids[i]];
              const auto& intF_shapeI = fe_values.intS_shapeI[f](i);

              const auto psi_bndry = EvaluateBoundaryCondition(bndry_id, node, groupset);

              for (size_t n = 0; n < num_gs_angles; ++n)
              {
                const auto& omega = quadrature->GetOmega(n);
                const auto mu = omega.Dot(face.normal);
                if (mu < 0.0)
                {
                  const auto& wt = quadrature->GetWeight(n);
                  const auto weight = -mu * wt * intF_shapeI;
                  const auto dof_map = discretization.MapDOFLocal(cell, i, uk_man, n, 0);

                  for (unsigned int gsg = 0; gsg < num_gs_groups; ++gsg)
                    local_response +=
                      weight * psi_dagger[gs][dof_map + gsg] * psi_bndry[num_gs_groups * n + gsg];
                } // if outgoing
              }
            } // for face node fi
          }
          ++f;
        } // for face
      } // for cell
      ++gs;
    } // for groupset
  } // if boundary sources

  // Point sources
  for (const auto& point_source : point_sources_)
    for (const auto& subscriber : point_source->GetSubscribers())
    {
      const auto& cell = grid->local_cells[subscriber.cell_local_id];
      const auto& transport_view = transport_views[cell.local_id];

      const auto src = point_source->GetStrength(0.0, num_groups);
      const auto& vol_wt = subscriber.volume_weight;

      const auto num_cell_nodes = transport_view.GetNumNodes();
      for (int i = 0; i < num_cell_nodes; ++i)
      {
        const auto dof_map = transport_view.MapDOF(i, 0, 0);
        const auto& shape_val = subscriber.shape_values(i);
        for (unsigned int g = 0; g < num_groups; ++g)
          local_response += vol_wt * shape_val * src[g] * phi_dagger[dof_map + g];
      } // for node i
    } // for subscriber

  // Volumetric sources
  for (const auto& volumetric_source : volumetric_sources_)
    for (const std::uint32_t local_id : volumetric_source->GetSubscribers())
    {
      const auto& cell = grid->local_cells[local_id];
      const auto& transport_view = transport_views[cell.local_id];
      const auto& fe_values = unit_cell_matrices[cell.local_id];
      const auto& nodes = discretization.GetCellNodeLocations(cell);

      const auto num_cell_nodes = transport_view.GetNumNodes();
      for (int i = 0; i < num_cell_nodes; ++i)
      {
        const auto& V_i = fe_values.intV_shapeI(i);
        const auto dof_map = transport_view.MapDOF(i, 0, 0);
        const auto& vals = volumetric_source->Evaluate(cell, nodes[i], num_groups, 0.0);
        for (unsigned int g = 0; g < num_groups; ++g)
          local_response += vals[g] * phi_dagger[dof_map + g] * V_i;
      }
    }

  double global_response = 0.0;
  mpi_comm.all_reduce(local_response, global_response, mpi::op::sum<double>());
  return global_response;
}

std::vector<double>
ResponseEvaluator::EvaluateBoundaryCondition(const uint64_t boundary_id,
                                             const Vector3& node,
                                             const LBSGroupset& groupset,
                                             const double /*unused*/) const
{
  const auto num_gs_angles = groupset.quadrature->GetNumAngles();
  const auto num_gs_groups = groupset.GetNumGroups();
  const auto first_group = groupset.first_group;

  std::vector<double> psi;
  const auto& bc = boundary_sources_.at(boundary_id);
  if (bc.type == LBSBoundaryType::ISOTROPIC)
  {
    for (size_t n = 0; n < num_gs_angles; ++n)
      for (unsigned int gsg = 0; gsg < num_gs_groups; ++gsg)
        psi.emplace_back(bc.isotropic_mg_source[first_group + gsg]);
    return psi;
  }
  OpenSnLogicalError("Unexpected behavior. Unsupported boundary condition encountered.");
}

} // namespace opensn

1	// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
2	// SPDX-License-Identifier: MIT
3
4	#include "modules/linear_boltzmann_solvers/response_evaluator/response_evaluator.h"
5	#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h"
6	#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/io/discrete_ordinates_problem_io.h"
7	#include "modules/linear_boltzmann_solvers/lbs_problem/point_source/point_source.h"
8	#include "modules/linear_boltzmann_solvers/lbs_problem/volumetric_source/volumetric_source.h"
9	#include "modules/linear_boltzmann_solvers/lbs_problem/io/lbs_problem_io.h"
10	#include "framework/mesh/mesh_continuum/mesh_continuum.h"
11	#include "framework/logging/log.h"
12	#include "framework/object_factory.h"
13	#include "framework/runtime.h"
14	#include "mpicpp-lite/mpicpp-lite.h"
15
16	namespace mpi = mpicpp_lite;
17
18	namespace opensn
19	{
20
21	OpenSnRegisterObjectInNamespace(lbs, ResponseEvaluator);
22
23	InputParameters
24	ResponseEvaluator::GetInputParameters()	32✔
25	{
26	InputParameters params;	32✔
27	params.SetGeneralDescription(	64✔
28	"A utility class for evaluating responses using precomputed adjoint solutions "
29	"and arbitrary forward sources.");
30
31	params.AddRequiredParameter<std::shared_ptr<Problem>>("problem",	64✔
32	"A handle to an existing LBS problem.");
33	params.AddOptionalParameterBlock(	64✔
34	"options", ParameterBlock(), "The specification of adjoint buffers and forward to use.");	64✔
35	params.LinkParameterToBlock("options", "response::OptionsBlock");	64✔
36
37	return params;	32✔
38	}	×
39
40	std::shared_ptr<ResponseEvaluator>
41	ResponseEvaluator::Create(const ParameterBlock& params)	32✔
42	{
43	auto& factory = opensn::ObjectFactory::GetInstance();	32✔
44	return factory.Create<ResponseEvaluator>("lbs::ResponseEvaluator", params);	64✔
45	}
46
47	ResponseEvaluator::ResponseEvaluator(const InputParameters& params)	32✔
48	: do_problem_(params.GetSharedPtrParam<Problem, DiscreteOrdinatesProblem>("problem"))	64✔
49	{
50	if (params.IsParameterValid("options"))	32✔
51	{
52	auto options = GetOptionsBlock();	×
53	options.AssignParameters(params.GetParam("options"));	×
54	SetOptions(options);	×
55	}	×
56	}	32✔
57
58	InputParameters
59	ResponseEvaluator::GetOptionsBlock()	32✔
60	{
61	InputParameters params;	32✔
62	params.SetGeneralDescription("A block of options for the response evaluator for adding adjoint "	64✔
63	"buffers and defining forward sources.");
64
65	params.AddOptionalParameterArray(	64✔
66	"buffers", {}, "An array of tables containing adjoint buffer specifications.");
67	params.LinkParameterToBlock("buffers", "response::BufferOptionsBlock");	64✔
68
69	params.AddOptionalParameter("clear_sources", false, "A flag to clear existing sources.");	64✔
70	params.AddOptionalParameterBlock(	64✔
71	"sources", ParameterBlock(), "An array of tables containing source specification information.");	64✔
72	params.LinkParameterToBlock("sources", "response::SourceOptionsBlock");	64✔
73
74	return params;	32✔
75	}	×
76
77	void
78	ResponseEvaluator::SetOptions(const InputParameters& params)	32✔
79	{
80	if (params.IsParameterValid("buffers"))	32✔
81	{
82	const auto& user_buffer_params = params.GetParam("buffers");	32✔
83	user_buffer_params.RequireBlockTypeIs(ParameterBlockType::ARRAY);	32✔
84	for (int p = 0; p < user_buffer_params.GetNumParameters(); ++p)	64✔
85	{
86	auto buffer_params = GetBufferOptionsBlock();	32✔
87	buffer_params.AssignParameters(user_buffer_params.GetParam(p));	32✔
88	SetBufferOptions(buffer_params);	32✔
89	}	32✔
90	}
91
92	if (params.IsParameterValid("clear_sources") and params.GetParamValue<bool>("clear_sources"))	32✔
93	{
94	material_sources_.clear();	×
95	point_sources_.clear();	×
96	volumetric_sources_.clear();	×
97	boundary_sources_.clear();	×
98	}
99
100	if (params.IsParameterValid("sources"))	32✔
101	{
102	auto source_params = GetSourceOptionsBlock();	32✔
103	source_params.AssignParameters(params.GetParam("sources"));	32✔
104	SetSourceOptions(source_params);	32✔
105	}	32✔
106	}	32✔
107
108	InputParameters
109	ResponseEvaluator::GetBufferOptionsBlock()	64✔
110	{
111	InputParameters params;	64✔
112	params.SetGeneralDescription("Options for adding adjoint buffers to the response evaluator.");	128✔
113
114	params.AddRequiredParameter<std::string>(	128✔
115	"name",
116	"A name given to the buffer to identify it when querying the response evaluation routine.");
117	params.AddRequiredParameterBlock(	128✔
118	"file_prefixes",
119	"A table containing file prefixes for flux moments and angular flux binary files. "
120	"These are keyed by \"flux_moments\" and \"angular_fluxes\", respectively.");
121
122	return params;	64✔
123	}	×
124
125	void
126	ResponseEvaluator::SetBufferOptions(const InputParameters& input)	32✔
127	{
128	auto params = opensn::ResponseEvaluator::GetBufferOptionsBlock();	32✔
129	params.AssignParameters(input);	32✔
130
131	const auto name = params.GetParamValue<std::string>("name");	32✔
132	OpenSnInvalidArgumentIf(adjoint_buffers_.count(name) > 0,	32✔
133	"An adjoint buffer with name " + name + " already exists.");
134
135	const auto prefixes = params.GetParam("file_prefixes");	32✔
136
137	std::vector<double> phi;	32✔
138	if (prefixes.Has("flux_moments"))	32✔
139	LBSSolverIO::ReadFluxMoments(	64✔
140	*do_problem_, prefixes.GetParamValue<std::string>("flux_moments"), false, phi);	64✔
141
142	std::vector<std::vector<double>> psi;	32✔
143	if (prefixes.Has("angular_fluxes"))	32✔
144	DiscreteOrdinatesProblemIO::ReadAngularFluxes(	×
145	*do_problem_, prefixes.GetParamValue<std::string>("angular_fluxes"), psi);	×
146
147	adjoint_buffers_[name] = {phi, psi};	32✔
148	log.Log0Verbose1() << "Adjoint buffer " << name << " added to the stack.";	64✔
149	}	32✔
150
151	InputParameters
152	ResponseEvaluator::GetSourceOptionsBlock()	64✔
153	{
154	InputParameters params;	64✔
155	params.SetGeneralDescription("A table of various forward source specifications.");	128✔
156
157	params.AddOptionalParameterArray(	128✔
158	"material", {}, "An array of tables containing material source specifications.");
159	params.LinkParameterToBlock("material", "response::MaterialSourceOptionsBlock");	128✔
160
161	params.AddOptionalParameterArray<std::shared_ptr<PointSource>>(	192✔
162	"point", {}, "An array of tables containing point source handles.");
163
164	params.AddOptionalParameterArray<std::shared_ptr<VolumetricSource>>(	128✔
165	"volumetric", {}, "An array of tables containing volumetric source handles.");
166
167	params.AddOptionalParameterArray(	128✔
168	"boundary", {}, "An array of tables containing boundary source specifications.");
169	params.LinkParameterToBlock("boundary", "response::BoundarySourceOptionsBlock");	128✔
170
171	return params;	64✔
172	}	×
173
174	InputParameters
175	ResponseEvaluator::GetBoundarySourceOptionsBlock()	×
176	{
177	InputParameters params;	×
178
179	params.SetGeneralDescription("Boundary source specification for response evaluations.");	×
180	params.SetClassName("Response Boundary Source");	×
181	params.AddRequiredParameter<std::string>("name",	×
182	"Boundary name that identifies the specific boundary");
183	params.AddRequiredParameter<std::string>("type", "Boundary type specification.");	×
184	params.AddOptionalParameterArray<double>("group_strength",	×
185	{},
186	"Required only if \"type\" is \"isotropic\". An array "
187	"of isotropic strength per group");
188	params.ConstrainParameterRange("type", AllowableRangeList::New({"isotropic"}));	×
189
190	return params;	×
191	}	×
192
193	void
194	ResponseEvaluator::SetSourceOptions(const InputParameters& input)	32✔
195	{
196	auto params = ResponseEvaluator::GetSourceOptionsBlock();	32✔
197	params.AssignParameters(input);	32✔
198
199	params.RequireBlockTypeIs(ParameterBlockType::BLOCK);	32✔
200
201	// Add material sources
202	if (params.Has("material"))	32✔
203	{
204	const auto& user_msrc_params = params.GetParam("material");	32✔
205	for (int p = 0; p < user_msrc_params.GetNumParameters(); ++p)	56✔
206	{
207	auto msrc_params = GetMaterialSourceOptionsBlock();	24✔
208	msrc_params.AssignParameters(user_msrc_params.GetParam(p));	24✔
209	SetMaterialSourceOptions(msrc_params);	24✔
210	}	24✔
211	}
212
213	// Add point sources
214	if (params.Has("point"))	32✔
215	{
216	const auto& user_psrc_params = params.GetParam("point");	32✔
217	for (int p = 0; p < user_psrc_params.GetNumParameters(); ++p)	40✔
218	{
219	point_sources_.push_back(	8✔
220	user_psrc_params.GetParam(p).GetValue<std::shared_ptr<PointSource>>());	8✔
221	point_sources_.back()->Initialize(*do_problem_);	8✔
222	}
223	}
224
225	// Add volumetric sources
226	if (params.Has("volumetric"))	32✔
227	{
228	const auto& user_dsrc_params = params.GetParam("volumetric");	32✔
229	for (int p = 0; p < user_dsrc_params.GetNumParameters(); ++p)	32✔
230	{
231	volumetric_sources_.push_back(	×
232	user_dsrc_params.GetParam(p).GetValue<std::shared_ptr<VolumetricSource>>());	×
233	volumetric_sources_.back()->Initialize(*do_problem_);	×
234	}
235	}
236
237	// Add boundary sources
238	if (params.Has("boundary"))	32✔
239	{
240	const auto& user_bsrc_params = params.GetParam("boundary");	32✔
241	for (int p = 0; p < user_bsrc_params.GetNumParameters(); ++p)	32✔
242	{
243	auto bsrc_params = ResponseEvaluator::GetBoundarySourceOptionsBlock();	×
244	bsrc_params.AssignParameters(user_bsrc_params.GetParam(p));	×
245	SetBoundarySourceOptions(bsrc_params);	×
246	}	×
247	}
248	}	32✔
249
250	InputParameters
251	ResponseEvaluator::GetMaterialSourceOptionsBlock()	24✔
252	{
253	InputParameters params;	24✔
254	params.SetGeneralDescription(	48✔
255	"Options for adding material-based forward sources to the response evaluator.");
256
257	params.AddRequiredParameter<unsigned int>("block_id", "The block id the source belongs to.");	48✔
258	params.AddRequiredParameterArray("strength", "The group-wise material source strength.");	48✔
259
260	return params;	24✔
261	}	×
262
263	void
264	ResponseEvaluator::SetMaterialSourceOptions(const InputParameters& params)	24✔
265	{
266	const auto blkid = params.GetParamValue<unsigned int>("block_id");	24✔
267	OpenSnInvalidArgumentIf(material_sources_.count(blkid) > 0,	24✔
268	"A material source for block id " + std::to_string(blkid) +
269	" already exists.");
270
271	const auto values = params.GetParamVectorValue<double>("strength");	24✔
272	OpenSnInvalidArgumentIf(values.size() != do_problem_->GetNumGroups(),	24✔
273	"The number of material source values and groups "
274	"in the underlying solver do not match. "
275	"Expected " +
276	std::to_string(do_problem_->GetNumGroups()) + " but got " +
277	std::to_string(values.size()) + ".");
278
279	material_sources_[blkid] = values;	24✔
280	log.Log0Verbose1() << "Material source for block id " << blkid << " added to the stack.";	48✔
281	}	24✔
282
283	void
284	ResponseEvaluator::SetBoundarySourceOptions(const InputParameters& params)	×
285	{
286	const auto bndry_name = params.GetParamValue<std::string>("name");	×
287	if (params.IsParameterValid("function"))	×
288	throw std::runtime_error("Boundary '" + bndry_name +	×
289	"' in ResponseEvaluator does not support \"function\".");	×
290	const auto bndry_type = params.GetParamValue<std::string>("type");	×
291
292	auto grid = do_problem_->GetGrid();	×
293	const auto bnd_name_map = grid->GetBoundaryNameMap();	×
294	OpenSnInvalidArgumentIf(not bnd_name_map.count(bndry_name),	×
295	"Boundary \"" + bndry_name + "\" does not exist.");
296	const auto bid = bnd_name_map.at(bndry_name);	×
297	if (bndry_type == "isotropic")	×
298	{
299	OpenSnInvalidArgumentIf(not params.Has("group_strength"),	×
300	"Parameter \"group_strength\" is required for "
301	"boundaries of type \"isotropic\".");
302	params.RequireParameterBlockTypeIs("group_strength", ParameterBlockType::ARRAY);	×
303	const auto values = params.GetParamVectorValue<double>("group_strength");	×
304	OpenSnInvalidArgumentIf(values.size() != do_problem_->GetNumGroups(),	×
305	"The number of boundary source values does not match the "
306	"number of groups. Expected " +
307	std::to_string(do_problem_->GetNumGroups()) + " but got " +
308	std::to_string(values.size()) + ".");
309
310	boundary_sources_[bid] = {LBSBoundaryType::ISOTROPIC, values};	×
311	}	×
312	else
313	log.Log0Warning() << "Unsupported boundary type. Skipping the entry.";	×
314	}	×
315
316	void
317	ResponseEvaluator::ClearForwardSources()	×
318	{
319	material_sources_.clear();	×
320	point_sources_.clear();	×
321	volumetric_sources_.clear();	×
322	boundary_sources_.clear();	×
323	}	×
324
325	void
326	ResponseEvaluator::AddResponseBuffers(const InputParameters& params)	×
327	{
328	params.RequireBlockTypeIs(ParameterBlockType::ARRAY);	×
329	for (size_t p = 0; p < params.GetNumParameters(); ++p)	×
330	{
331	auto spec = ResponseEvaluator::GetBufferOptionsBlock();	×
332	spec.AssignParameters(params.GetParam(p));	×
333	SetBufferOptions(spec);	×
334	}	×
335	}	×
336
337	void
338	ResponseEvaluator::AddResponseSources(const InputParameters& params)	×
339	{
340	auto spec = ResponseEvaluator::GetSourceOptionsBlock();	×
341	spec.AssignParameters(params);	×
342	SetSourceOptions(spec);	×
343	}	×
344
345	double
346	ResponseEvaluator::EvaluateResponse(const std::string& buffer) const	32✔
347	{
348	const auto& buffer_data = adjoint_buffers_.at(buffer);	32✔
349	const auto& phi_dagger = buffer_data.first;	32✔
350	const auto& psi_dagger = buffer_data.second;	32✔
351
352	OpenSnLogicalErrorIf(not material_sources_.empty() and phi_dagger.empty(),	32✔
353	"If material sources are present, adjoint flux moments "
354	"must be available for response evaluation.");
355	OpenSnLogicalErrorIf(not point_sources_.empty() and phi_dagger.empty(),	32✔
356	"If point sources are set, adjoint flux moments "
357	"must be available for response evaluation.");
358	OpenSnLogicalErrorIf(not volumetric_sources_.empty() and phi_dagger.empty(),	32✔
359	"if volumetric sources are set, adjoint flux moments "
360	"must be available for response evaluation.");
361	OpenSnLogicalErrorIf(not boundary_sources_.empty() and psi_dagger.empty(),	32✔
362	"If boundary sources are set, adjoint angular fluxes "
363	"must be available for response evaluation.");
364
365	const auto& grid = do_problem_->GetGrid();	32✔
366	const auto& discretization = do_problem_->GetSpatialDiscretization();	32✔
367	const auto& transport_views = do_problem_->GetCellTransportViews();	32✔
368	const auto& unit_cell_matrices = do_problem_->GetUnitCellMatrices();	32✔
369	const auto num_groups = do_problem_->GetNumGroups();	32✔
370
371	double local_response = 0.0;	32✔
372
373	// Material sources
374	if (not material_sources_.empty())	32✔
375	{
376	for (const auto& cell : grid->local_cells)	16,898✔
377	{
378	const auto& cell_mapping = discretization.GetCellMapping(cell);	16,874✔
379	const auto& transport_view = transport_views[cell.local_id];	16,874✔
380	const auto& fe_values = unit_cell_matrices[cell.local_id];	16,874✔
381	const auto num_cell_nodes = cell_mapping.GetNumNodes();	16,874✔
382
383	if (material_sources_.count(cell.block_id) > 0)	33,442✔
384	{
385	const auto& src = material_sources_.at(cell.block_id);	306✔
386	for (size_t i = 0; i < num_cell_nodes; ++i)	1,050✔
387	{
388	const auto dof_map = transport_view.MapDOF(i, 0, 0);	744✔
389	const auto& V_i = fe_values.intV_shapeI(i);	744✔
390	for (unsigned int g = 0; g < num_groups; ++g)	16,176✔
391	local_response += src[g] * phi_dagger[dof_map + g] * V_i;	15,432✔
392	}
393	}
394	} // for cell
395	} // if material sources
396
397	// Boundary sources
398	if (not boundary_sources_.empty())	32✔
399	{
400	size_t gs = 0;	×
401	for (const auto& groupset : do_problem_->GetGroupsets())	×
402	{
403	const auto& uk_man = groupset.psi_uk_man_;	×
404	const auto& quadrature = groupset.quadrature;	×
NEW 405	const auto num_gs_angles = quadrature->GetNumAngles();	×
406	const auto& num_gs_groups = groupset.GetNumGroups();	×
407
408	for (const auto& cell : grid->local_cells)	×
409	{
410	const auto& cell_mapping = discretization.GetCellMapping(cell);	×
411	const auto& fe_values = unit_cell_matrices[cell.local_id];	×
412
413	size_t f = 0;	×
414	for (const auto& face : cell.faces)	×
415	{
416	if (not face.has_neighbor and boundary_sources_.count(face.neighbor_id) > 0)	×
417	{
418	const auto bndry_id = face.neighbor_id;	×
419	const auto num_face_nodes = cell_mapping.GetNumFaceNodes(f);	×
420	for (size_t fi = 0; fi < num_face_nodes; ++fi)	×
421	{
422	const auto i = cell_mapping.MapFaceNode(f, fi);	×
423	const auto& node = grid->vertices[cell.vertex_ids[i]];	×
424	const auto& intF_shapeI = fe_values.intS_shapeI[f](i);	×
425
426	const auto psi_bndry = EvaluateBoundaryCondition(bndry_id, node, groupset);	×
427
428	for (size_t n = 0; n < num_gs_angles; ++n)	×
429	{
NEW 430	const auto& omega = quadrature->GetOmega(n);	×
431	const auto mu = omega.Dot(face.normal);	×
432	if (mu < 0.0)	×
433	{
NEW 434	const auto& wt = quadrature->GetWeight(n);	×
435	const auto weight = -mu * wt * intF_shapeI;	×
436	const auto dof_map = discretization.MapDOFLocal(cell, i, uk_man, n, 0);	×
437
438	for (unsigned int gsg = 0; gsg < num_gs_groups; ++gsg)	×
439	local_response +=	×
440	weight * psi_dagger[gs][dof_map + gsg] * psi_bndry[num_gs_groups * n + gsg];	×
441	} // if outgoing
442	}
443	} // for face node fi	×
444	}
445	++f;	×
446	} // for face
447	} // for cell
448	++gs;	×
449	} // for groupset
450	} // if boundary sources
451
452	// Point sources
453	for (const auto& point_source : point_sources_)	40✔
454	for (const auto& subscriber : point_source->GetSubscribers())	10✔
455	{
456	const auto& cell = grid->local_cells[subscriber.cell_local_id];	2✔
457	const auto& transport_view = transport_views[cell.local_id];	2✔
458
459	const auto src = point_source->GetStrength(0.0, num_groups);	2✔
460	const auto& vol_wt = subscriber.volume_weight;	2✔
461
462	const auto num_cell_nodes = transport_view.GetNumNodes();	2✔
463	for (int i = 0; i < num_cell_nodes; ++i)	10✔
464	{
465	const auto dof_map = transport_view.MapDOF(i, 0, 0);	8✔
466	const auto& shape_val = subscriber.shape_values(i);	8✔
467	for (unsigned int g = 0; g < num_groups; ++g)	52✔
468	local_response += vol_wt * shape_val * src[g] * phi_dagger[dof_map + g];	44✔
469	} // for node i
470	} // for subscriber	2✔
471
472	// Volumetric sources
473	for (const auto& volumetric_source : volumetric_sources_)	32✔
474	for (const std::uint32_t local_id : volumetric_source->GetSubscribers())	×
475	{
476	const auto& cell = grid->local_cells[local_id];	×
477	const auto& transport_view = transport_views[cell.local_id];	×
478	const auto& fe_values = unit_cell_matrices[cell.local_id];	×
479	const auto& nodes = discretization.GetCellNodeLocations(cell);	×
480
481	const auto num_cell_nodes = transport_view.GetNumNodes();	×
482	for (int i = 0; i < num_cell_nodes; ++i)	×
483	{
484	const auto& V_i = fe_values.intV_shapeI(i);	×
485	const auto dof_map = transport_view.MapDOF(i, 0, 0);	×
486	const auto& vals = volumetric_source->Evaluate(cell, nodes[i], num_groups, 0.0);	×
487	for (unsigned int g = 0; g < num_groups; ++g)	×
488	local_response += vals[g] * phi_dagger[dof_map + g] * V_i;	×
489	}	×
490	}
491
492	double global_response = 0.0;	32✔
493	mpi_comm.all_reduce(local_response, global_response, mpi::op::sum<double>());	32✔
494	return global_response;	32✔
495	}	32✔
496
497	std::vector<double>
498	ResponseEvaluator::EvaluateBoundaryCondition(const uint64_t boundary_id,	×
499	const Vector3& node,
500	const LBSGroupset& groupset,
501	const double /unused/) const
502	{
NEW 503	const auto num_gs_angles = groupset.quadrature->GetNumAngles();	×
504	const auto num_gs_groups = groupset.GetNumGroups();	×
505	const auto first_group = groupset.first_group;	×
506
507	std::vector<double> psi;	×
508	const auto& bc = boundary_sources_.at(boundary_id);	×
509	if (bc.type == LBSBoundaryType::ISOTROPIC)	×
510	{
511	for (size_t n = 0; n < num_gs_angles; ++n)	×
512	for (unsigned int gsg = 0; gsg < num_gs_groups; ++gsg)	×
513	psi.emplace_back(bc.isotropic_mg_source[first_group + gsg]);	×
514	return psi;	×
515	}
516	OpenSnLogicalError("Unexpected behavior. Unsupported boundary condition encountered.");	×
517	}	×
518
519	} // namespace opensn

Open-Sn / opensn / 25654184358

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous