Skip to content

Commit

Permalink
Merge Enrico Degregori's PICO improvements into dev
Browse files Browse the repository at this point in the history
  • Loading branch information
@ckhroulev
ckhroulev committed Jun 9, 2021
2 parents e67ba13 + 7b505c5 commit 58dfa19
Show file tree
Hide file tree
Showing 4 changed files with 73 additions and 37 deletions.
73 changes: 43 additions & 30 deletions src/coupler/ocean/Pico.cc
View file
Expand Up @@ -39,7 +39,6 @@
#include "pism/util/iceModelVec.hh"
#include "pism/util/Time.hh"
#include "pism/geometry/Geometry.hh"

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

#include "Pico.hh"
Expand Down Expand Up @@ -328,8 +327,11 @@ 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) {

std::vector<int> count(m_n_basins, 0);
// additional vectors to allreduce efficiently with IntelMPI
std::vector<int> countr(m_n_basins, 0);
std::vector<double> salinityr(m_n_basins);
std::vector<double> temperaturer(m_n_basins);

temperature.resize(m_n_basins);
salinity.resize(m_n_basins);
Expand Down Expand Up @@ -359,11 +361,16 @@ 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.
for (int basin_id = 0; basin_id < m_n_basins; basin_id++) {
GlobalSum(m_grid->com, count.data(), countr.data(), m_n_basins);
GlobalSum(m_grid->com, salinity.data(), salinityr.data(), m_n_basins);
GlobalSum(m_grid->com, temperature.data(), temperaturer.data(), m_n_basins);

// copy values
count = countr;
salinity = salinityr;
temperature = temperaturer;

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]);
for (int basin_id = 0; basin_id < m_n_basins; 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 Down Expand Up @@ -399,11 +406,14 @@ 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) {

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<int> n_shelf_cells_per_basin(m_n_shelves*m_n_basins,0);
std::vector<int> n_shelf_cells(m_n_shelves, 0);
// additional vectors to allreduce efficiently with IntelMPI
std::vector<int> n_shelf_cells_per_basinr(m_n_shelves*m_n_basins,0);
std::vector<int> n_shelf_cellsr(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
Expand All @@ -412,16 +422,15 @@ void Pico::set_ocean_input_fields(const PicoPhysics &physics, const IceModelVec2
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*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_cellsr.data(), m_n_shelves);
GlobalSum(m_grid->com, n_shelf_cells_per_basin.data(), n_shelf_cells_per_basinr.data(), m_n_shelves*m_n_basins);
// copy data
n_shelf_cells = n_shelf_cellsr;
n_shelf_cells_per_basin = n_shelf_cells_per_basinr;
}

// now set potential temperature and salinity box 0:
Expand All @@ -441,8 +450,8 @@ 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*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 Down Expand Up @@ -749,7 +758,7 @@ 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);
// fill results with zeros
result.resize(m_n_shelves);
for (int s = 0; s < m_n_shelves; ++s) {
Expand All @@ -767,15 +776,17 @@ void Pico::compute_box_average(int box_id,
result[shelf_id] += field(i, j);
}
}

// compute the global sum and average
for (int s = 0; s < m_n_shelves; ++s) {
auto n_cells = GlobalSum(m_grid->com, n_cells_per_box[s]);
std::vector<int> n_cells(m_n_shelves);
GlobalSum(m_grid->com, n_cells_per_box.data(), n_cells.data(), m_n_shelves);
GlobalSum(m_grid->com, result.data(), result1.data(), m_n_shelves);

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

if (n_cells > 0) {
result[s] /= (double)n_cells;
for (int s = 0; s < m_n_shelves; ++s) {
if (n_cells[s] > 0) {
result[s] /= (double)n_cells[s];
}
}
}
Expand All @@ -793,7 +804,6 @@ void Pico::compute_box_area(int box_id,
const IceModelVec2Int &box_mask,
std::vector<double> &result) {
result.resize(m_n_shelves);

IceModelVec::AccessList list{ &shelf_mask, &box_mask };

auto cell_area = m_grid->cell_area();
Expand All @@ -807,10 +817,13 @@ 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]);

// compute GlobalSum from index 1 to index m_n_shelves-1
std::vector<double> result1(m_n_shelves);
GlobalSum(m_grid->com, &result[1], &result1[1], m_n_shelves-1);
// copy data
for (int i = 1; i < m_n_shelves; i++) {
result[i] = result1[i];
}
}

Expand Down
21 changes: 14 additions & 7 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,7 +18,6 @@
*/

#include <algorithm> // max_element

#include "PicoGeometry.hh"
#include "pism/util/connected_components.hh"
#include "pism/util/IceModelVec2CellType.hh"
Expand Down Expand Up @@ -139,6 +138,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 +161,13 @@ static void relabel(RelabelingType type,
loop.failed();
}
loop.check();
GlobalSum(grid->com, area.data(), area1.data(), area.size());

// copy data
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 +627,9 @@ 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);

for (Points p(*m_grid); p; p.next()) {
const int i = p.i(), j = p.j();
Expand All @@ -647,10 +653,11 @@ 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);
// copy data
GL_distance_max = GL_distance_max1;
CF_distance_max = CF_distance_max1;

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

Expand Down
12 changes: 12 additions & 0 deletions src/util/pism_utilities.cc
View file
Expand Up @@ -22,6 +22,7 @@
#include <cstdarg> // va_list, va_start(), va_end()
#include <sstream> // istringstream, ostringstream
#include <cstdio> // vsnprintf
#include <cassert> // assert

#include <mpi.h> // MPI_Get_library_version
#include <fftw3.h> // fftw_version
Expand Down Expand Up @@ -128,10 +129,17 @@ bool member(const std::string &string, const std::set<std::string> &set) {
}

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

void GlobalReduce(MPI_Comm comm, int *local, int *result, int count, MPI_Op op) {
assert(local != result);
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 +152,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
4 changes: 4 additions & 0 deletions src/util/pism_utilities.hh
View file
Expand Up @@ -107,12 +107,16 @@ double vector_max(const std::vector<double> &input);
// parallel
void GlobalReduce(MPI_Comm comm, double *local, double *result, int count, MPI_Op op);

void GlobalReduce(MPI_Comm comm, int *local, int *result, int count, MPI_Op op);

void GlobalMin(MPI_Comm comm, double *local, double *result, int count);

void GlobalMax(MPI_Comm comm, double *local, double *result, int count);

void GlobalSum(MPI_Comm comm, double *local, double *result, int count);

void GlobalSum(MPI_Comm comm, int *local, int *result, int count);

double GlobalMin(MPI_Comm comm, double local);

double GlobalMax(MPI_Comm comm, double local);
Expand Down

0 comments on commit 58dfa19

Please sign in to comment.