Skip to content

Fix multithreading reproducibility #98

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

Open
lrobion wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Fix multithreading reproducibility #98

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
81a619b
Fix race conditions
lrobion May 10, 2026
77a61fa
Disable parallel reductions
lrobion May 12, 2026
44a44b8
Cleanup compiler warnings
lrobion May 14, 2026
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
40 changes: 31 additions & 9 deletions Code.v05-00/src/AIM/Aerosol.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1128,8 +1128,10 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) {
Vector_2D moment(Ny, Vector_1D(Nx, 0.0E+00));
const double FACTOR = 3.0 / double(4.0 * PI);

/*
#pragma omp parallel for default(shared) private(iNx, jNy, iBin) \
schedule(dynamic, 1) if (!PARALLEL_CASES)
*/
for (iBin = 0; iBin < nBin; iBin++)
{
for (jNy = 0; jNy < Ny; jNy++)
Expand Down Expand Up @@ -1189,10 +1191,12 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) {
UInt iNx = 0;
UInt jNy = 0;

/*
#pragma omp parallel for default(shared) private(iNx, jNy) \
reduction(+ \
reduction(+ \
: totalnumber_sum) \
schedule(dynamic, 1) if (!PARALLEL_CASES)
schedule(dynamic, 1) if (!PARALLEL_CASES)
*/
for (jNy = 0; jNy < Ny; jNy++)
{
for (iNx = 0; iNx < Nx; iNx++)
Expand All @@ -1213,8 +1217,10 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) {
UInt iNx = 0;
UInt iBin = 0;

/*
#pragma omp parallel for default(shared) private(iNx, jNy, iBin) \
schedule(dynamic, 1) if (!PARALLEL_CASES)
schedule(dynamic, 1) if (!PARALLEL_CASES)
*/
for (iBin = 0; iBin < nBin; iBin++)
{
for (jNy = 0; jNy < Ny; jNy++)
Expand Down Expand Up @@ -1312,10 +1318,12 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) {
UInt iNx = 0;
UInt jNy = 0;

/*
#pragma omp parallel for default(shared) private(iNx, jNy) \
reduction(+ \
: totalicemass_sum) \
: totalicemass_sum) \
schedule(dynamic, 1) if (!PARALLEL_CASES)
*/
for (jNy = 0; jNy < Ny; jNy++)
{
for (iNx = 0; iNx < Nx; iNx++)
Expand All @@ -1335,8 +1343,10 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) {
Vector_2D TVol = TotalVolume();
const double FACTOR = RHO_ICE * 1.0E+06;

/*
#pragma omp parallel for default(shared) private(iNx, jNy) \
schedule(dynamic, 1) if (!PARALLEL_CASES)
schedule(dynamic, 1) if (!PARALLEL_CASES)
*/
for (jNy = 0; jNy < Ny; jNy++)
{
for (iNx = 0; iNx < Nx; iNx++)
Expand Down Expand Up @@ -1436,8 +1446,10 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) {

Vector_1D PDF(nBin, 0.0E+00);

/*
#pragma omp parallel for default(shared) private(iNx, jNy, iBin) \
schedule(dynamic, 1) if (!PARALLEL_CASES)
schedule(dynamic, 1) if (!PARALLEL_CASES)
*/
for (iBin = 0; iBin < nBin; iBin++)
{
for (jNy = 0; jNy < Ny; jNy++)
Expand Down Expand Up @@ -1466,8 +1478,10 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) {

Vector_1D PDF(nBin, 0.0E+00);

/*
#pragma omp parallel for default(shared) private(iNx, jNy, iBin) \
schedule(dynamic, 1) if (!PARALLEL_CASES)
schedule(dynamic, 1) if (!PARALLEL_CASES)
*/
for (iBin = 0; iBin < nBin; iBin++)
{
for (jNy = 0; jNy < Ny; jNy++)
Expand Down Expand Up @@ -1619,10 +1633,12 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) {

double moment = 0.0E+00;

/*
#pragma omp parallel for default(shared) private(iBin) \
reduction(+ \
reduction(+ \
: moment) \
schedule(dynamic, 1) if (!PARALLEL_CASES)
*/
for (iBin = 0; iBin < nBin; iBin++)
{
moment += (log(bin_Edges[iBin + 1] / bin_Edges[iBin])) * pow(bin_Centers[iBin], n) * PDF[iBin];
Expand All @@ -1640,10 +1656,12 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) {
double moment = 0.0E+00;
const double FACTOR = 3.0 / double(4.0 * PI);

/*
#pragma omp parallel for default(shared) private(iBin) \
reduction(+ \
: moment) \
schedule(dynamic, 1) if (!PARALLEL_CASES)
*/
for (iBin = 0; iBin < nBin; iBin++)
moment += (log(bin_Edges[iBin + 1] / bin_Edges[iBin])) * pow(FACTOR * bin_VCenters[iBin][jNy][iNx], n / double(3.0)) * pdf[iBin][jNy][iNx];

Expand Down Expand Up @@ -1714,8 +1732,10 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) {

Vector_1D PDF(nBin, 0.0E+00);

/*
#pragma omp parallel for default(shared) private(iNx, jNy, iBin) \
schedule(dynamic, 1) if (!PARALLEL_CASES)
*/
for (iBin = 0; iBin < nBin; iBin++)
{
for (jNy = 0; jNy < Ny; jNy++)
Expand Down Expand Up @@ -1758,8 +1778,10 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) {
UInt jNy = 0;
UInt iBin = 0;

/*
#pragma omp parallel for default(shared) private(iNx, jNy, iBin) \
schedule(dynamic, 1) if (!PARALLEL_CASES)
schedule(dynamic, 1) if (!PARALLEL_CASES)
*/
for (jNy = 0; jNy < Ny; jNy++)
{
for (iNx = 0; iNx < Nx; iNx++)
Expand Down
8 changes: 4 additions & 4 deletions Code.v05-00/src/Core/LAGRIDPlumeModel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -187,10 +187,10 @@ SimStatus LAGRIDPlumeModel::runFullModel() {
// dp/dz = -rho*g = -(n/V)Mg
Vector_3D& pdfRef = iceAerosol_.getPDF();
auto pressureEdges = met_.PressEdges();
double localND;
// double localND;
#pragma omp parallel for
for (std::size_t j=0; j<yCoords_.size(); j++){
localND = (pressureEdges[j] - pressureEdges[j+1])/(yEdges_[j+1] - yEdges_[j]);
double localND = (pressureEdges[j] - pressureEdges[j+1])/(yEdges_[j+1] - yEdges_[j]);
for (std::size_t i=0; i<xCoords_.size(); i++){
Contrail_[j][i] = Contrail_[j][i] / localND; // parts per trillion
H2O_[j][i] = H2O_[j][i] / localND; // parts per trillion
Expand Down Expand Up @@ -763,10 +763,10 @@ void LAGRIDPlumeModel::remapAllVars(double remapTimestep, const std::vector<std:
// included it). This uses the hydrostatic assumption:
// dp/dz = -rho*g = -(n/V)Mg
auto pressureEdges = met_.PressEdges();
double localND;
// double localND;
#pragma omp parallel for
for (std::size_t j=0; j<yCoords_.size(); j++){
localND = (pressureEdges[j] - pressureEdges[j+1])/(yEdges_[j+1] - yEdges_[j]);
double localND = (pressureEdges[j] - pressureEdges[j+1])/(yEdges_[j+1] - yEdges_[j]);
for (std::size_t i=0; i<xCoords_.size(); i++){
Contrail_[j][i] = Contrail_[j][i] * localND;
H2O_[j][i] = H2O_[j][i] * localND;
Expand Down
6 changes: 4 additions & 2 deletions Code.v05-00/src/Core/Meteorology.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -506,9 +506,11 @@ Vector_1D Meteorology::interpMetTimeseriesData(double simTime_h, const Vector_2D
}

void Meteorology::updateTempPerturb() {
/*
#pragma omp parallel for\
if(!PARALLEL_CASES) \
default(shared)
if(!PARALLEL_CASES) \
default(shared)
*/
for (int j = 0; j < ny_; j++){
for(int i = 0; i < nx_; i++){
double epsilon1 = fRand(-1.0, 1.0);
Expand Down