Skip to content

Split solution vectors per groupset #666

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
andrsd wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

Split solution vectors per groupset #666

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions framework/math/math.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,13 +63,29 @@ Scale(std::vector<double>& x, const double& val)
xi *= val;
}

void
Scale(std::vector<std::vector<double>>& xs, const double& val)
{
for (auto& x : xs)
for (double& xi : x)
xi *= val;
}

void
Set(std::vector<double>& x, const double& val)
{
for (double& xi : x)
xi = val;
}

void
Set(std::vector<std::vector<double>>& xs, const double& val)
{
for (auto& x : xs)
for (double& xi : x)
xi = val;
}

double
Dot(const std::vector<double>& x, const std::vector<double>& y)
{
Expand Down
4 changes: 4 additions & 0 deletions framework/math/math.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,9 +95,13 @@ void PrintVector(const std::vector<double>& x);
/// Scales a vector in place by constant.
void Scale(std::vector<double>& x, const double& val);

void Scale(std::vector<std::vector<double>>& xs, const double& val);

/// Sets a constant value to a vector.
void Set(std::vector<double>& x, const double& val);

void Set(std::vector<std::vector<double>>& xss, const double& val);

/// Multiplies the vector with a constant and returns result.
std::vector<double> Mult(const std::vector<double>& x, const double& val);

Expand Down
2 changes: 2 additions & 0 deletions framework/utils/hdf_utils.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,10 @@
#pragma once

#include "hdf5.h"
#include <stdexcept>
#include <vector>
#include <string>
#include <cassert>

namespace opensn
{
Expand Down
4 changes: 2 additions & 2 deletions ..._solvers/discrete_ordinates_curvilinear_problem/discrete_ordinates_curvilinear_problem.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -365,9 +365,9 @@ DiscreteOrdinatesCurvilinearProblem::SetSweepChunk(LBSGroupset& groupset)
secondary_unit_cell_matrices_,
cell_transport_views_,
densities_local_,
phi_new_local_,
phi_new_local_[groupset.id],
psi_new_local_[groupset.id],
q_moments_local_,
q_moments_local_[groupset.id],
groupset,
block_id_to_xs_map_,
num_moments_,
Expand Down
10 changes: 5 additions & 5 deletions ...ltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/aah_sweep_chunk_rz.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ AAHSweepChunkRZ::AAHSweepChunkRZ(const std::shared_ptr<MeshContinuum> grid,
const SpatialDiscretization& discretization_primary,
const std::vector<UnitCellMatrices>& unit_cell_matrices,
const std::vector<UnitCellMatrices>& secondary_unit_cell_matrices,
std::vector<CellLBSView>& cell_transport_views,
std::vector<std::vector<CellLBSView>>& cell_transport_views,
const std::vector<double>& densities,
std::vector<double>& destination_phi,
std::vector<double>& destination_psi,
Expand Down Expand Up @@ -98,7 +98,7 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set)
auto cell_local_id = spls[spls_index];
auto& cell = grid_->local_cells[cell_local_id];
auto& cell_mapping = discretization_.GetCellMapping(cell);
auto& cell_transport_view = cell_transport_views_[cell_local_id];
auto& cell_transport_view = cell_transport_views_[groupset_.id][cell_local_id];
auto cell_num_faces = cell.faces.size();
auto cell_num_nodes = cell_mapping.GetNumNodes();

Expand Down Expand Up @@ -247,7 +247,7 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set)
double temp_src = 0.0;
for (int m = 0; m < num_moments_; ++m)
{
const size_t ir = cell_transport_view.MapDOF(i, m, static_cast<int>(gs_gi + gsg));
const size_t ir = cell_transport_view.MapDOF(i, m, static_cast<int>(gsg));
temp_src += m2d_op[m][direction_num] * source_moments_[ir];
}
source[i] = temp_src;
Expand Down Expand Up @@ -278,7 +278,7 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set)
const double wn_d2m = d2m_op[m][direction_num];
for (int i = 0; i < cell_num_nodes; ++i)
{
const size_t ir = cell_transport_view.MapDOF(i, m, gs_gi);
const size_t ir = cell_transport_view.MapDOF(i, m, 0);
for (int gsg = 0; gsg < gs_size; ++gsg)
destination_phi_[ir + gsg] += wn_d2m * b[gsg](i);
}
Expand Down Expand Up @@ -327,7 +327,7 @@ AAHSweepChunkRZ::Sweep(AngleSet& angle_set)
{
for (int gsg = 0; gsg < gs_size; ++gsg)
cell_transport_view.AddOutflow(
f, gs_gi + gsg, wt * face_mu_values[f] * b[gsg](i) * IntF_shapeI(i));
f, gsg, wt * face_mu_values[f] * b[gsg](i) * IntF_shapeI(i));
}

double* psi = nullptr;
Expand Down
2 changes: 1 addition & 1 deletion ...oltzmann_solvers/discrete_ordinates_curvilinear_problem/sweep_chunks/aah_sweep_chunk_rz.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ class AAHSweepChunkRZ : public SweepChunk
const SpatialDiscretization& discretization_primary,
const std::vector<UnitCellMatrices>& unit_cell_matrices,
const std::vector<UnitCellMatrices>& secondary_unit_cell_matrices,
std::vector<CellLBSView>& cell_transport_views,
std::vector<std::vector<CellLBSView>>& cell_transport_views,
const std::vector<double>& densities,
std::vector<double>& destination_phi,
std::vector<double>& destination_psi,
Expand Down
137 changes: 69 additions & 68 deletions modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -182,13 +182,13 @@ DiscreteOrdinatesProblem::ReorientAdjointSolution()

for (const auto& groupset : groupsets_)
{
int gs = groupset.id;
int gsi = groupset.id;

// Moment map for flux moments
const auto& moment_map = groupset.quadrature->GetMomentToHarmonicsIndexMap();

// Angular flux info
auto& psi = psi_new_local_[gs];
auto& psi = psi_new_local_[gsi];
const auto& uk_man = groupset.psi_uk_man_;

// Build reversed angle mapping
Expand Down Expand Up @@ -228,18 +228,16 @@ DiscreteOrdinatesProblem::ReorientAdjointSolution()

OpenSnLogicalErrorIf(not found,
"Opposing angle for " + omegas[idir].PrintStr() + " in groupset " +
std::to_string(gs) + " not found.");
std::to_string(gsi) + " not found.");

} // for angle m
} // if saving angular flux

const auto num_gs_groups = groupset.groups.size();
const auto gsg_i = groupset.groups.front().id;
const auto gsg_f = groupset.groups.back().id;

for (const auto& cell : grid_->local_cells)
{
const auto& transport_view = cell_transport_views_[cell.local_id];
const auto& transport_view = cell_transport_views_[gsi][cell.local_id];
for (int i = 0; i < transport_view.GetNumNodes(); ++i)
{
// Reorient flux moments
Expand All @@ -255,10 +253,10 @@ DiscreteOrdinatesProblem::ReorientAdjointSolution()
const auto& ell = moment_map[imom].ell;
const auto dof_map = transport_view.MapDOF(i, imom, 0);

for (int g = gsg_i; g <= gsg_f; ++g)
for (int g = 0; g < num_gs_groups; ++g)
{
phi_new_local_[dof_map + g] *= std::pow(-1.0, ell);
phi_old_local_[dof_map + g] *= std::pow(-1.0, ell);
phi_new_local_[gsi][dof_map + g] *= std::pow(-1.0, ell);
phi_old_local_[gsi][dof_map + g] *= std::pow(-1.0, ell);
} // for group g
} // for moment m

Expand All @@ -277,8 +275,7 @@ DiscreteOrdinatesProblem::ReorientAdjointSolution()
}
} // for node i
} // for cell

} // for groupset
} // for groupset
}

void
Expand All @@ -288,8 +285,8 @@ DiscreteOrdinatesProblem::ZeroOutflowBalanceVars(LBSGroupset& groupset)

for (const auto& cell : grid_->local_cells)
for (int f = 0; f < cell.faces.size(); ++f)
for (auto& group : groupset.groups)
cell_transport_views_[cell.local_id].ZeroOutflow(f, group.id);
for (std::size_t g = 0; g < groupset.groups.size(); ++g)
cell_transport_views_[groupset.id][cell.local_id].ZeroOutflow(f, g);
}

void
Expand All @@ -302,52 +299,55 @@ DiscreteOrdinatesProblem::ComputeBalance()
// Get material source
// This is done using the SetSource routine because it allows a lot of flexibility.
auto mat_src = phi_new_local_;
mat_src.assign(mat_src.size(), 0.0);
Set(mat_src, 0.);
for (auto& groupset : groupsets_)
{
q_moments_local_.assign(q_moments_local_.size(), 0.0);
Set(q_moments_local_, 0.);
active_set_source_function_(groupset,
q_moments_local_,
phi_new_local_,
APPLY_FIXED_SOURCES | APPLY_AGS_FISSION_SOURCES |
APPLY_WGS_FISSION_SOURCES);
LBSVecOps::GSScopedCopyPrimarySTLvectors(*this, groupset, q_moments_local_, mat_src);
LBSVecOps::GSScopedCopyPrimarySTLvectors(
*this, groupset, q_moments_local_[groupset.id], mat_src[groupset.id]);
}

// Compute absorption, material-source and in-flow
double local_out_flow = 0.0;
double local_in_flow = 0.0;
double local_absorption = 0.0;
double local_production = 0.0;
for (const auto& cell : grid_->local_cells)
for (auto& groupset : groupsets_)
{
const auto& cell_mapping = discretization_->GetCellMapping(cell);
const auto& transport_view = cell_transport_views_[cell.local_id];
const auto& fe_intgrl_values = unit_cell_matrices_[cell.local_id];
const size_t num_nodes = transport_view.GetNumNodes();
const auto& IntV_shapeI = fe_intgrl_values.intV_shapeI;
const auto& IntS_shapeI = fe_intgrl_values.intS_shapeI;

// Inflow: This is essentially an integration over all faces, all angles, and all groups. For
// non-reflective boundaries, only the cosines that are negative are added to the inflow
// integral. For reflective boundaries, it is expected that, upon convergence, inflow = outflow
// (within numerical tolerances set by the user).
for (int f = 0; f < cell.faces.size(); ++f)
const auto gsi = groupset.id;
const auto first_group = groupset.groups.front().id;
for (const auto& cell : grid_->local_cells)
{
const auto& face = cell.faces[f];

if (not face.has_neighbor) // Boundary face
const auto& cell_mapping = discretization_->GetCellMapping(cell);
const auto& transport_view = cell_transport_views_[gsi][cell.local_id];
const auto& fe_intgrl_values = unit_cell_matrices_[cell.local_id];
const size_t num_nodes = cell_mapping.GetNumNodes();
const auto& IntV_shapeI = fe_intgrl_values.intV_shapeI;
const auto& IntS_shapeI = fe_intgrl_values.intS_shapeI;

// Inflow: This is essentially an integration over all faces, all angles, and all groups. For
// non-reflective boundaries, only the cosines that are negative are added to the inflow
// integral. For reflective boundaries, it is expected that, upon convergence, inflow =
// outflow (within numerical tolerances set by the user).
for (int f = 0; f < cell.faces.size(); ++f)
{
const auto& bndry = sweep_boundaries_[face.neighbor_id];
const auto& face = cell.faces[f];

if (bndry->IsReflecting())
{
for (int g = 0; g < num_groups_; ++g)
local_in_flow += transport_view.GetOutflow(f, g);
}
else
if (not face.has_neighbor) // Boundary face
{
for (const auto& groupset : groupsets_)
const auto& bndry = sweep_boundaries_[face.neighbor_id];

if (bndry->IsReflecting())
{
for (int g = 0; g < groupset.groups.size(); ++g)
local_in_flow += transport_view.GetOutflow(f, g);
}
else
{
for (int n = 0; n < groupset.quadrature->omegas.size(); ++n)
{
Expand All @@ -371,32 +371,33 @@ DiscreteOrdinatesProblem::ComputeBalance()
} // for fi
} // if mu < 0
} // for n
} // for groupset
} // if reflecting boundary
} // if boundary
} // for f

// Outflow: The group-wise outflow was determined during a solve so we just accumulate it here.
for (int f = 0; f < cell.faces.size(); ++f)
for (int g = 0; g < num_groups_; ++g)
local_out_flow += transport_view.GetOutflow(f, g);

// Absorption and sources
const auto& xs = transport_view.GetXS();
const auto& sigma_a = xs.GetSigmaAbsorption();
for (int i = 0; i < num_nodes; ++i)
{
for (int g = 0; g < num_groups_; ++g)
} // if reflecting boundary
} // if boundary
} // for f

// Outflow: The group-wise outflow was determined during a solve so we just accumulate it
// here.
for (int f = 0; f < cell.faces.size(); ++f)
for (int g = 0; g < groupset.groups.size(); ++g)
local_out_flow += transport_view.GetOutflow(f, g);

// Absorption and sources
const auto& xs = transport_view.GetXS();
const auto& sigma_a = xs.GetSigmaAbsorption();
for (int i = 0; i < num_nodes; ++i)
{
size_t imap = transport_view.MapDOF(i, 0, g);
double phi_0g = phi_new_local_[imap];
double q_0g = mat_src[imap];
for (int g = 0; g < groupset.groups.size(); ++g)
{
size_t imap = transport_view.MapDOF(i, 0, g);
double phi_0g = phi_new_local_[gsi][imap];
double q_0g = mat_src[gsi][imap];

local_absorption += sigma_a[g] * phi_0g * IntV_shapeI(i);
local_production += q_0g * IntV_shapeI(i);
} // for g
} // for i
} // for cell
local_absorption += sigma_a[first_group + g] * phi_0g * IntV_shapeI(i);
local_production += q_0g * IntV_shapeI(i);
} // for g
} // for i
} // for cell
} // groupsets

// Compute local balance
double local_balance = local_production + local_in_flow - local_absorption - local_out_flow;
Expand Down Expand Up @@ -1076,9 +1077,9 @@ DiscreteOrdinatesProblem::SetSweepChunk(LBSGroupset& groupset)
unit_cell_matrices_,
cell_transport_views_,
densities_local_,
phi_new_local_,
phi_new_local_[groupset.id],
psi_new_local_[groupset.id],
q_moments_local_,
q_moments_local_[groupset.id],
groupset,
block_id_to_xs_map_,
num_moments_,
Expand All @@ -1088,14 +1089,14 @@ DiscreteOrdinatesProblem::SetSweepChunk(LBSGroupset& groupset)
}
else if (sweep_type_ == "CBC")
{
auto sweep_chunk = std::make_shared<CBCSweepChunk>(phi_new_local_,
auto sweep_chunk = std::make_shared<CBCSweepChunk>(phi_new_local_[groupset.id],
psi_new_local_[groupset.id],
grid_,
*discretization_,
unit_cell_matrices_,
cell_transport_views_,
densities_local_,
q_moments_local_,
q_moments_local_[groupset.id],
groupset,
block_id_to_xs_map_,
num_moments_,
Expand Down
Loading
Loading