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

MHD reproducibility #648

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
+105 −307
Open

MHD reproducibility #648

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
6ddfb52
Minor fixes
tomidakn Feb 12, 2025
30faabf
Improved reproducibility of MHD simulations
tomidakn Feb 14, 2025
cdd766e
trailing whitespace
felker Feb 14, 2025
83b393a
style fix
tomidakn Feb 14, 2025
4622a99
Merge branch 'mhd_reproducibility' of https://github.com/PrincetonUni…
tomidakn Feb 14, 2025
1802abe
Fixed a bug in shearing box.
tomidakn Feb 14, 2025
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
1 change: 1 addition & 0 deletions src/bvals/fc/bvals_fc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,6 +162,7 @@ FaceCenteredBoundaryVariable::FaceCenteredBoundaryVariable(MeshBlock *pmb,
#endif
}
} // end "if is a shearing boundary"
ClearEMFShearing(shear_var_emf_[upper]);
} // end loop over inner, outer shearing boundaries
} // end shearing box component of ctor
}
Expand Down
1 change: 0 additions & 1 deletion src/bvals/fc/bvals_fc.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -189,7 +189,6 @@ class FaceCenteredBoundaryVariable : public BoundaryVariable {
void SetFluxBoundaryFromPolar(Real **buf_list, int num_bufs, bool is_north);

void ClearCoarseFluxBoundary();
void AverageFluxBoundary(); //!> average flux from fine and equal lvls
void PolarFluxBoundarySingleAzimuthalBlock();
void CountFineEdges(); //!> called in SetupPersistentMPI()

Expand Down
33 changes: 9 additions & 24 deletions src/bvals/fc/bvals_shear_emf.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
//========================================================================================

// C headers
#include <float.h>
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
#include <float.h>
#include <float.h>

see below, use #include <limits> (or at a minimum, cfloat)


// C++ headers
#include <algorithm> // min
Expand Down Expand Up @@ -100,28 +101,7 @@ void FaceCenteredBoundaryVariable::SendEMFShearingBoxBoundaryCorrection() {
int offset[2]{0, 3};
for (int upper=0; upper<2; upper++) {
if (pbval_->is_shear[upper]) {
// step 1. -- average edges of shboxvar_fc_flx_
// average e3 for x1x2 edge
for (int k=pmb->ks; k<=pmb->ke; k++) {
for (int j=pmb->js; j<=pmb->je+1; j+=pmb->block_size.nx2)
shear_var_emf_[upper].x3e(k,j) *= 0.5;
}
if(pmb->block_size.nx3 > 1) {
// average e2 for x1x3 edge
if (pbval_->block_bcs[BoundaryFace::inner_x3] == BoundaryFlag::block
|| pbval_->block_bcs[BoundaryFace::inner_x3] == BoundaryFlag::periodic) {
for (int j=pmb->js; j<=pmb->je; j++)
shear_var_emf_[upper].x2e(pmb->ks,j) *= 0.5;
}
if (pbval_->block_bcs[BoundaryFace::outer_x3] == BoundaryFlag::block
|| pbval_->block_bcs[BoundaryFace::outer_x3] == BoundaryFlag::periodic) {
for (int j=pmb->js; j<=pmb->je; j++)
shear_var_emf_[upper].x2e(pmb->ke+1,j) *= 0.5;
}
}

// step 2. -- load sendbuf; memcpy to recvbuf if on same rank, post
// MPI_Isend otherwise
// load sendbuf; memcpy to recvbuf if on same rank, post MPI_Isend otherwise
for (int n=0; n<3; n++) {
SimpleNeighborBlock& snb = pbval_->sb_flux_data_[upper].send_neighbor[n];
if (snb.rank != -1) {
Expand Down Expand Up @@ -298,7 +278,12 @@ void FaceCenteredBoundaryVariable::SetEMFShearingBoxBoundaryCorrection() {
void FaceCenteredBoundaryVariable::ClearEMFShearing(EdgeField &work) {
AthenaArray<Real> &e2 = work.x2e;
AthenaArray<Real> &e3 = work.x3e;
e2.ZeroClear();
e3.ZeroClear();
constexpr Real m = (sizeof(Real) == sizeof(float)) ? -FLT_MAX : -DBL_MAX;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the C++ way of doing this is much better:

constexpr Real m = -std::numeric_limits<Real>::max();

I noticed that the C-style approach somehow made it through to orbital advection and FC flux correction files in the past:

❯ grep -nri "flt_max" ./
./orbital_advection/default_orbital_velocity.cpp:11:#include <cfloat>     // FLT_MAX, FLT_MIN
./orbital_advection/orbital_remapping.cpp:11:#include <cfloat>     // FLT_MAX
./orbital_advection/set_orbital_advection.cpp:11:#include <cfloat>     // FLT_MAX
./orbital_advection/calculate_orbital_advection.cpp:11:#include <cfloat>     // FLT_MAX
./orbital_advection/orbital_advection.cpp:11:#include <cfloat>     // FLT_MAX, FLT_MIN
./orbital_advection/orbital_advection.cpp:435:    vK_max = -(FLT_MAX);
./orbital_advection/orbital_advection.cpp:436:    vK_min = (FLT_MAX);
./orbital_advection/orbital_advection.cpp:512:    min_dt = (FLT_MAX);
./bvals/fc/flux_correction_fc.cpp:1224:  constexpr Real m = (sizeof(Real) == sizeof(float)) ? -FLT_MAX : -DBL_MAX;

would you mind changing them?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please change them. Regarding OA, I just opened another PR about it, so probably it is better to fix them there.

int s = e2.GetSize();
for (int i = 0; i < s; ++i)
e2(i) = m;
s = e3.GetSize();
for (int i = 0; i < s; ++i)
e3(i) = m;
return;
}
Loading