Skip to content
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

WIP: Update cell density based on cell-avg'ed temperature #135

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

WIP: Update cell density based on cell-avg'ed temperature #135

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: 1 addition & 15 deletions include/enrico/coupled_driver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,12 +49,7 @@ class CoupledDriver {
//! Update the temperature for the neutronics solver
//!
//! \param relax Apply relaxation to temperature before updating neutronics solver
void update_temperature(bool relax);

//! Update the density for the neutronics solver
//!
//! \param relax Apply relaxation to density before updating neutronics solver
void update_density(bool relax);
void update_temperature_and_density(bool relax);

//! Check convergence of the coupled solve for the current Picard iteration.
bool is_converged();
Expand Down Expand Up @@ -170,15 +165,6 @@ class CoupledDriver {

xt::xtensor<double, 1> temperatures_prev_; //!< Previous Picard iteration temperature

//! Current Picard iteration density; this density is the density
//! computed by the thermal-hydraulic solver, and data mappings may result in
//! a different density actually used in the neutronics solver. For example,
//! the entries in this xtensor may be averaged over neutronics cells to give
//! the density used by the neutronics solver.
xt::xtensor<double, 1> densities_;

xt::xtensor<double, 1> densities_prev_; //!< Previous Picard iteration density

//! Current Picard iteration heat source; this heat source is the heat source
//! computed by the neutronics solver, and data mappings may result in a different
//! heat source actually used in the heat solver. For example, the entries in this
Expand Down
127 changes: 8 additions & 119 deletions src/coupled_driver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#include "enrico/comm_split.h"
#include "enrico/driver.h"
#include "enrico/error.h"
#include "iapws/iapws.h"

#ifdef USE_NEK5000
#include "enrico/nek5000_driver.h"
Expand Down Expand Up @@ -94,18 +95,6 @@ CoupledDriver::CoupledDriver(MPI_Comm comm, pugi::xml_node node)
}
}

if (coup_node.child("density_ic")) {
std::string s = coup_node.child_value("density_ic");

if (s == "neutronics") {
density_ic_ = Initial::neutronics;
} else if (s == "heat_fluids") {
density_ic_ = Initial::heat;
} else {
throw std::runtime_error{"Invalid value for <density_ic>"};
}
}

Expects(power_ > 0);
Expects(max_timesteps_ >= 0);
Expects(max_picard_iter_ >= 0);
Expand Down Expand Up @@ -186,7 +175,6 @@ CoupledDriver::CoupledDriver(MPI_Comm comm, pugi::xml_node node)
init_cell_fluid_mask();

init_temperatures();
init_densities();
init_heat_source();
}

Expand Down Expand Up @@ -232,8 +220,7 @@ void CoupledDriver::execute()
// At this point, there is always a previous iterate of temperature and density
// (as assured by the initial conditions set in init_temperature and init_density)
// so we always apply underrelaxation here.
update_temperature(true);
update_density(true);
update_temperature_and_density(true);

if (is_converged()) {
std::string msg = "converged at i_picard = " + std::to_string(i_picard_);
Expand Down Expand Up @@ -332,9 +319,9 @@ void CoupledDriver::update_heat_source(bool relax)
}
}

void CoupledDriver::update_temperature(bool relax)
void CoupledDriver::update_temperature_and_density(bool relax)
{
comm_.message("Updating temperature");
comm_.message("Updating temperature and density");

auto& neutronics = this->get_neutronics_driver();
auto& heat = this->get_heat_driver();
Expand Down Expand Up @@ -386,59 +373,11 @@ void CoupledDriver::update_temperature(bool relax)
// Set temperature for cell instance
CellHandle cell = kv.first;
neutronics.set_temperature(cell, average_temp);
}
}
}

void CoupledDriver::update_density(bool relax)
{

auto& neutronics = this->get_neutronics_driver();
auto& heat = this->get_heat_driver();

// *********************************************************************
// Gather densities on heat root and apply underrelaxation
// *********************************************************************

if (relax && comm_.rank == heat_root_) {
std::copy(densities_.begin(), densities_.end(), densities_prev_.begin());
}

densities_ = heat.density();

if (relax && comm_.rank == heat_root_) {
if (alpha_rho_ == ROBBINS_MONRO) {
int n = i_picard_ + 1;
densities_ = densities_ / n + (1. - 1. / n) * densities_prev_;
} else {
densities_ = alpha_rho_ * densities_ + (1.0 - alpha_rho_) * densities_prev_;
}
}

// ************************************************************************
// Send underrelaxed density to all neutron ranks and set cell temperatures
// ************************************************************************

comm_.send_and_recv(densities_, neutronics_root_, heat_root_);
neutronics.comm_.broadcast(densities_);

if (neutronics.active()) {
// For each neutronics cell in a fluid cell, volume average the densities
// and set in driver
for (const auto& kv : cell_to_elems_) {
CellHandle cell = kv.first;
// If cell is in fluid, set density for cell instance
if (cell_fluid_mask_[cell] == 1) {
// Get volume-averaged density for this cell instance
double average_density = 0.0;
double total_vol = 0.0;
for (int e : kv.second) {
average_density += densities_[e] * elem_volumes_[e];
total_vol += elem_volumes_[e];
}
// Set density for cell instance
average_density /= total_vol;
Ensures(average_density > 0.0);
neutronics.set_density(cell, average_density);
double density = 1e-3 / iapws::nu1(heat.pressure_bc_, average_temp);
neutronics.set_density(cell, density);
}
}
}
Expand Down Expand Up @@ -521,7 +460,7 @@ void CoupledDriver::init_temperatures()
// temperatures received from the heat solver.
// * We do not want to apply underrelaxation here (and at this point,
// there is no previous iterate of temperature, anyway).
update_temperature(false);
update_temperature_and_density(false);
}

// In both cases, only temperatures_ was set, so we explicitly set temperatures_prev_
Expand Down Expand Up @@ -561,56 +500,6 @@ void CoupledDriver::init_volumes()
}
}

void CoupledDriver::init_densities()
{
comm_.message("Initializing densities");

const auto& neutronics = this->get_neutronics_driver();
const auto& heat = this->get_heat_driver();

if (comm_.rank == heat_root_) {
densities_.resize({static_cast<unsigned long>(n_global_elem_)});
densities_prev_.resize({static_cast<unsigned long>(n_global_elem_)});
} else if (comm_.rank == neutronics_root_) {
densities_.resize({static_cast<unsigned long>(n_global_elem_)});
}

if (density_ic_ == Initial::neutronics) {
if (comm_.rank == neutronics_root_) {
// Loop over the neutronics cells, then loop over the TH elements
// corresponding to that cell and assign the cell density to the correct
// index in the densities_ array. This mapping assumes that each TH
// element is fully contained within a neutronics cell, i.e., TH elements
// are not split between multiple neutronics cells.
for (const auto& kv : cell_to_elems_) {
CellHandle cell = kv.first;
const auto& global_elems = kv.second;

for (int elem : global_elems) {
if (cell_fluid_mask_[cell] == 1) {
double rho = neutronics.get_density(cell);
densities_[elem] = rho;
} else {
densities_[elem] = 0.0;
}
}
}
}
comm_.send_and_recv(densities_, heat_root_, neutronics_root_);
} else if (density_ic_ == Initial::heat) {
// * This sets densities_ on the the coupling_root, based on the
// densities received from the heat solver.
// * We do not want to apply underrelaxation here (and at this point,
// there is no previous iterate of density, anyway).
update_density(false);
}

// In both cases, we need to explicitly set densities_prev_
if (comm_.rank == heat_root_) {
std::copy(densities_.begin(), densities_.end(), densities_prev_.begin());
}
}

void CoupledDriver::init_elem_fluid_mask()
{
comm_.message("Initializing element fluid mask");
Expand Down