Skip to content

Commit

Permalink
PICO optimized communication
Browse files Browse the repository at this point in the history
  • Loading branch information
@ckhroulev
k202136 authored and ckhroulev committed Jun 9, 2021
1 parent e67ba13 commit 5718be7
Show file tree
Hide file tree
Showing 4 changed files with 129 additions and 37 deletions.
117 changes: 88 additions & 29 deletions src/coupler/ocean/Pico.cc
View file
Expand Up @@ -39,7 +39,7 @@
#include "pism/util/iceModelVec.hh"
#include "pism/util/Time.hh"
#include "pism/geometry/Geometry.hh"

#include "pism/util/Profiling.hh"
#include "pism/coupler/util/options.hh"

#include "Pico.hh"
Expand Down Expand Up @@ -328,8 +328,12 @@ void Pico::compute_ocean_input_per_basin(const PicoPhysics &physics, const IceMo
const IceModelVec2Int &continental_shelf_mask,
const IceModelVec2S &salinity_ocean, const IceModelVec2S &theta_ocean,
std::vector<double> &temperature, std::vector<double> &salinity) {

const Profiling &profiling = m_grid->ctx()->profiling();
profiling.begin("compute_ocean_input_per_basin");
std::vector<int> count(m_n_basins, 0);
std::vector<int> count1(m_n_basins, 0);
std::vector<double> salinity1(m_n_basins);
std::vector<double> temperature1(m_n_basins);

temperature.resize(m_n_basins);
salinity.resize(m_n_basins);
Expand Down Expand Up @@ -359,11 +363,18 @@ void Pico::compute_ocean_input_per_basin(const PicoPhysics &physics, const IceMo
// ocean_contshelf_mask values intersect with the basin, count is zero. In such case,
// use dummy temperature and salinity. This could happen, for example, if the ice shelf
// front advances beyond the continental shelf break.
GlobalSum(m_grid->com, count.data(), count1.data(), m_n_basins);
GlobalSum(m_grid->com, salinity.data(), salinity1.data(), m_n_basins);
GlobalSum(m_grid->com, temperature.data(), temperature1.data(), m_n_basins);
count = count1;
salinity = salinity1;
temperature = temperature1;

for (int basin_id = 0; basin_id < m_n_basins; basin_id++) {

count[basin_id] = GlobalSum(m_grid->com, count[basin_id]);
salinity[basin_id] = GlobalSum(m_grid->com, salinity[basin_id]);
temperature[basin_id] = GlobalSum(m_grid->com, temperature[basin_id]);
// count[basin_id] = GlobalSum(m_grid->com, count[basin_id]);
// salinity[basin_id] = GlobalSum(m_grid->com, salinity[basin_id]);
// temperature[basin_id] = GlobalSum(m_grid->com, temperature[basin_id]);

// if basin is not dummy basin 0 or there are no ocean cells in this basin to take the mean over.
if (basin_id > 0 && count[basin_id] == 0) {
Expand All @@ -385,6 +396,7 @@ void Pico::compute_ocean_input_per_basin(const PicoPhysics &physics, const IceMo
m_log->message(5, " %d: temp =%.3f, salinity=%.3f\n", basin_id, temperature[basin_id], salinity[basin_id]);
}
}
profiling.end("compute_ocean_input_per_basin");
}

//! Set ocean ocean input from box 0 as boundary condition for box 1.
Expand All @@ -399,29 +411,43 @@ void Pico::set_ocean_input_fields(const PicoPhysics &physics, const IceModelVec2
const IceModelVec2Int &shelf_mask, const std::vector<double> basin_temperature,
const std::vector<double> basin_salinity, IceModelVec2S &Toc_box0,
IceModelVec2S &Soc_box0) {


const Profiling &profiling = m_grid->ctx()->profiling();
profiling.begin("set_ocean_input_fields");
IceModelVec::AccessList list{ &ice_thickness, &basin_mask, &Soc_box0, &Toc_box0, &mask, &shelf_mask };

std::vector<std::vector<int> > n_shelf_cells_per_basin(m_n_shelves, std::vector<int>(m_n_basins, 0));
// std::vector<std::vector<int> > n_shelf_cells_per_basin(m_n_shelves, std::vector<int>(m_n_basins, 0));
std::vector<int> n_shelf_cells_per_basin(m_n_shelves*m_n_basins,0);
std::vector<int> n_shelf_cells(m_n_shelves, 0);
std::vector<int> n_shelf_cells_per_basin1(m_n_shelves*m_n_basins,0);
std::vector<int> n_shelf_cells1(m_n_shelves, 0);

// 1) count the number of cells in each shelf
// 2) count the number of cells in the intersection of each shelf with all the basins
{
int* ptr_shelf_cells;
int* ptr_shelf_cells_per_basin;

for (Points p(*m_grid); p; p.next()) {
const int i = p.i(), j = p.j();
int s = shelf_mask.as_int(i, j);
int b = basin_mask.as_int(i, j);
n_shelf_cells_per_basin[s][b]++;
// n_shelf_cells_per_basin[s][b]++;
n_shelf_cells_per_basin[s*m_n_basins+b]++;
n_shelf_cells[s]++;
}

for (int s = 0; s < m_n_shelves; s++) {
n_shelf_cells[s] = GlobalSum(m_grid->com, n_shelf_cells[s]);
for (int b = 0; b < m_n_basins; b++) {
n_shelf_cells_per_basin[s][b] = GlobalSum(m_grid->com, n_shelf_cells_per_basin[s][b]);
}
}

GlobalSum(m_grid->com, n_shelf_cells.data(), n_shelf_cells1.data(), m_n_shelves);
GlobalSum(m_grid->com, n_shelf_cells_per_basin.data(), n_shelf_cells_per_basin1.data(), m_n_shelves*m_n_basins);
n_shelf_cells = n_shelf_cells1;
n_shelf_cells_per_basin = n_shelf_cells_per_basin1;

// for (int s = 0; s < m_n_shelves; s++) {
// n_shelf_cells[s] = GlobalSum(m_grid->com, n_shelf_cells[s]);
// for (int b = 0; b < m_n_basins; b++) {
// n_shelf_cells_per_basin[s][b] = GlobalSum(m_grid->com, n_shelf_cells_per_basin[s][b]);
// }
// }
}

// now set potential temperature and salinity box 0:
Expand All @@ -441,8 +467,10 @@ void Pico::set_ocean_input_fields(const PicoPhysics &physics, const IceModelVec2

// weighted input depending on the number of shelf cells in each basin
for (int b = 1; b < m_n_basins; b++) { //Note: b=0 yields nan
Toc_box0(i, j) += basin_temperature[b] * n_shelf_cells_per_basin[s][b] / (double)n_shelf_cells[s];
Soc_box0(i, j) += basin_salinity[b] * n_shelf_cells_per_basin[s][b] / (double)n_shelf_cells[s];
// Toc_box0(i, j) += basin_temperature[b] * n_shelf_cells_per_basin[s][b] / (double)n_shelf_cells[s];
// Soc_box0(i, j) += basin_salinity[b] * n_shelf_cells_per_basin[s][b] / (double)n_shelf_cells[s];
Toc_box0(i, j) += basin_temperature[b] * n_shelf_cells_per_basin[s*m_n_basins+b] / (double)n_shelf_cells[s];
Soc_box0(i, j) += basin_salinity[b] * n_shelf_cells_per_basin[s*m_n_basins+b] / (double)n_shelf_cells[s];
}

double theta_pm = physics.theta_pm(Soc_box0(i, j), physics.pressure(ice_thickness(i, j)));
Expand All @@ -462,6 +490,7 @@ void Pico::set_ocean_input_fields(const PicoPhysics &physics, const IceModelVec2
"setting it to pressure melting temperature\n",
low_temperature_counter);
}
profiling.end("set_ocean_input_fields");
}

/*!
Expand Down Expand Up @@ -528,7 +557,8 @@ void Pico::process_box1(const PicoPhysics &physics,
IceModelVec2S &Toc,
IceModelVec2S &Soc,
IceModelVec2S &overturning) {

const Profiling &profiling = m_grid->ctx()->profiling();
profiling.begin("process_box1");
std::vector<double> box1_area(m_n_shelves);

compute_box_area(1, shelf_mask, box_mask, box1_area);
Expand Down Expand Up @@ -581,6 +611,7 @@ void Pico::process_box1(const PicoPhysics &physics,
"has been negative in %d cases!\n",
n_Toc_failures);
}
profiling.end("process_box1");
}

void Pico::process_other_boxes(const PicoPhysics &physics,
Expand All @@ -593,6 +624,8 @@ void Pico::process_other_boxes(const PicoPhysics &physics,
IceModelVec2S &Toc,
IceModelVec2S &Soc) {

const Profiling &profiling = m_grid->ctx()->profiling();
profiling.begin("process_other_boxes");
std::vector<double> overturning(m_n_shelves, 0.0);
std::vector<double> salinity(m_n_shelves, 0.0);
std::vector<double> temperature(m_n_shelves, 0.0);
Expand Down Expand Up @@ -667,6 +700,7 @@ void Pico::process_other_boxes(const PicoPhysics &physics,
}

} // loop over boxes
profiling.end("process_other_boxes");
}

/*!
Expand Down Expand Up @@ -749,7 +783,10 @@ void Pico::compute_box_average(int box_id,
IceModelVec::AccessList list{ &field, &shelf_mask, &box_mask };

std::vector<int> n_cells_per_box(m_n_shelves, 0);

std::vector<double> result1(m_n_shelves);
double* ptrres;
int* ptrlocal;
int* ptrresult;
// fill results with zeros
result.resize(m_n_shelves);
for (int s = 0; s < m_n_shelves; ++s) {
Expand All @@ -767,17 +804,30 @@ void Pico::compute_box_average(int box_id,
result[shelf_id] += field(i, j);
}
}

const Profiling &profiling = m_grid->ctx()->profiling();
// compute the global sum and average
profiling.begin("deg_test");
std::vector<int> n_cells(m_n_shelves);
ptrlocal = n_cells_per_box.data();
ptrresult = n_cells.data();
GlobalSum(m_grid->com, ptrlocal, ptrresult, m_n_shelves);
GlobalSum(m_grid->com, result.data(), result1.data(), m_n_shelves);
result = result1;
for (int s = 0; s < m_n_shelves; ++s) {
auto n_cells = GlobalSum(m_grid->com, n_cells_per_box[s]);
// auto n_cells = GlobalSum(m_grid->com, n_cells_per_box[s]);

result[s] = GlobalSum(m_grid->com, result[s]);
// result[s] = GlobalSum(m_grid->com, result[s]);

if (n_cells > 0) {
result[s] /= (double)n_cells;
if (n_cells[s] > 0) {
result[s] /= (double)n_cells[s];
}
// if (n_cells > 0) {
// result[s] /= (double)n_cells;
// }
}


profiling.end("deg_test");
}

/*!
Expand All @@ -793,9 +843,10 @@ void Pico::compute_box_area(int box_id,
const IceModelVec2Int &box_mask,
std::vector<double> &result) {
result.resize(m_n_shelves);

std::vector<double> result1(m_n_shelves);
IceModelVec::AccessList list{ &shelf_mask, &box_mask };

double* ptrres;
double* ptrres1;
auto cell_area = m_grid->cell_area();

for (Points p(*m_grid); p; p.next()) {
Expand All @@ -807,11 +858,19 @@ void Pico::compute_box_area(int box_id,
result[shelf_id] += cell_area;
}
}

// compute global sums
for (int s = 1; s < m_n_shelves; ++s) {
result[s] = GlobalSum(m_grid->com, result[s]);
ptrres = result.data(); // point to first value (index 0)
ptrres1 = result1.data();
ptrres++; // point to second value (index 1)
ptrres1++;
GlobalSum(m_grid->com, ptrres, ptrres1, m_n_shelves-1); // GlobalSum from index 1 to index m_n_shelves-1
for (int i = 1; i < m_n_shelves; i++) {
result[i] = result1[i];
}
// for (int s = 1; s < m_n_shelves; ++s) {
// result[s] = GlobalSum(m_grid->com, result[s]);
// }
}

} // end of namespace ocean
Expand Down
36 changes: 28 additions & 8 deletions src/coupler/ocean/PicoGeometry.cc
View file
@@ -1,4 +1,4 @@
/* Copyright (C) 2018, 2019, 2020 PISM Authors
/* Copyright (C) 2018, 2019, 2020, 2021 PISM Authors
*
* This file is part of PISM.
*
Expand All @@ -18,12 +18,13 @@
*/

#include <algorithm> // max_element

#include <numeric> // accumula6te vector
#include "PicoGeometry.hh"
#include "pism/util/connected_components.hh"
#include "pism/util/IceModelVec2CellType.hh"
#include "pism/util/pism_utilities.hh"
#include "pism/util/petscwrappers/Vec.hh"
#include "pism/util/Profiling.hh"

namespace pism {
namespace ocean {
Expand Down Expand Up @@ -77,6 +78,9 @@ const IceModelVec2Int &PicoGeometry::ice_rise_mask() const {
* to date.
*/
void PicoGeometry::update(const IceModelVec2S &bed_elevation, const IceModelVec2CellType &cell_type) {
const Profiling &profiling = m_grid->ctx()->profiling();
profiling.begin("ocean.update_geometry");
profiling.begin("ocean.update_geometry.1");
bool exclude_ice_rises = m_config->get_flag("ocean.pico.exclude_ice_rises");

int n_boxes = m_config->get_number("ocean.pico.number_of_boxes");
Expand All @@ -91,8 +95,10 @@ void PicoGeometry::update(const IceModelVec2S &bed_elevation, const IceModelVec2

compute_lakes(cell_type, m_lake_mask);
}
profiling.end("ocean.update_geometry.1");

// these two could be optimized by trying to reduce the number of times we update ghosts
profiling.begin("ocean.update_geometry.2");
{
m_ice_rises.update_ghosts();
m_ocean_mask.update_ghosts();
Expand All @@ -101,15 +107,19 @@ void PicoGeometry::update(const IceModelVec2S &bed_elevation, const IceModelVec2

compute_distances_cf(m_ocean_mask, m_ice_rises, exclude_ice_rises, m_distance_cf);
}
profiling.end("ocean.update_geometry.2");

// these two could be done at the same time
profiling.begin("ocean.update_geometry.3");
{
compute_ice_shelf_mask(m_ice_rises, m_lake_mask, m_ice_shelves);

compute_continental_shelf_mask(bed_elevation, m_ice_rises, continental_shelf_depth, m_continental_shelf);
}

compute_box_mask(m_distance_gl, m_distance_cf, m_ice_shelves, n_boxes, m_boxes);
profiling.end("ocean.update_geometry.3");
profiling.end("ocean.update_geometry");
}


Expand Down Expand Up @@ -139,6 +149,7 @@ static void relabel(RelabelingType type,
}

std::vector<double> area(max_index + 1, 0.0);
std::vector<double> area1(max_index + 1, 0.0);
{

ParallelSection loop(grid->com);
Expand All @@ -161,9 +172,10 @@ static void relabel(RelabelingType type,
loop.failed();
}
loop.check();

GlobalSum(grid->com, area.data(), area1.data(), area.size());
area = area1;
for (unsigned int k = 0; k < area.size(); ++k) {
area[k] = grid->cell_area() * GlobalSum(grid->com, area[k]);
area[k] = grid->cell_area() * area[k];
}
}

Expand Down Expand Up @@ -623,7 +635,11 @@ void PicoGeometry::compute_box_mask(const IceModelVec2Int &D_gl, const IceModelV
int n_shelves = shelf_mask.range().max + 1;

std::vector<double> GL_distance_max(n_shelves, 0.0);
std::vector<double> GL_distance_max1(n_shelves, 0.0);
std::vector<double> CF_distance_max(n_shelves, 0.0);
std::vector<double> CF_distance_max1(n_shelves, 0.0);
double* ptrGL;
double* ptrCF;

for (Points p(*m_grid); p; p.next()) {
const int i = p.i(), j = p.j();
Expand All @@ -647,10 +663,14 @@ void PicoGeometry::compute_box_mask(const IceModelVec2Int &D_gl, const IceModelV
}

// compute global maximums
for (int k = 0; k < n_shelves; ++k) {
GL_distance_max[k] = GlobalMax(m_grid->com, GL_distance_max[k]);
CF_distance_max[k] = GlobalMax(m_grid->com, CF_distance_max[k]);
}
GlobalMax(m_grid->com, GL_distance_max.data(), GL_distance_max1.data(), n_shelves);
GlobalMax(m_grid->com, CF_distance_max.data(), CF_distance_max1.data(), n_shelves);
GL_distance_max = GL_distance_max1;
CF_distance_max = CF_distance_max1;
// for (int k = 0; k < n_shelves; ++k) {
// GL_distance_max[k] = GlobalMax(m_grid->com, GL_distance_max[k]);
// CF_distance_max[k] = GlobalMax(m_grid->com, CF_distance_max[k]);
// }

double GL_distance_ref = *std::max_element(GL_distance_max.begin(), GL_distance_max.end());

Expand Down
9 changes: 9 additions & 0 deletions src/util/pism_utilities.cc
View file
Expand Up @@ -132,6 +132,11 @@ void GlobalReduce(MPI_Comm comm, double *local, double *result, int count, MPI_O
PISM_C_CHK(err, 0, "MPI_Allreduce");
}

void GlobalReduce(MPI_Comm comm, int *local, int *result, int count, MPI_Op op) {
int err = MPI_Allreduce(local, result, count, MPI_INT, op, comm);
PISM_C_CHK(err, 0, "MPI_Allreduce");
}

void GlobalMin(MPI_Comm comm, double *local, double *result, int count) {
GlobalReduce(comm, local, result, count, MPI_MIN);
}
Expand All @@ -144,6 +149,10 @@ void GlobalSum(MPI_Comm comm, double *local, double *result, int count) {
GlobalReduce(comm, local, result, count, MPI_SUM);
}

void GlobalSum(MPI_Comm comm, int *local, int *result, int count) {
GlobalReduce(comm, local, result, count, MPI_SUM);
}

unsigned int GlobalSum(MPI_Comm comm, unsigned int input) {
unsigned int result;
int err = MPI_Allreduce(&input, &result, 1, MPI_UNSIGNED, MPI_SUM, comm);
Expand Down

0 comments on commit 5718be7

Please sign in to comment.