Skip to content

Commit

Permalink
Fixed a bug in shearing box.
Browse files Browse the repository at this point in the history
  • Loading branch information
@tomidakn
tomidakn committed Feb 14, 2025
1 parent 4622a99 commit 1802abe
Show file tree
Hide file tree
Showing 3 changed files with 29 additions and 35 deletions.
1 change: 1 addition & 0 deletions src/bvals/fc/bvals_fc.cpp
View file
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]);

This comment has been minimized.

Sign in to view

Copy link
@changgoo

changgoo Feb 14, 2025

Contributor

oops, do you think this can introduce some bugs in the original calculations? I'm asking that I have some troubles in my simulations when shearing box turned on.

This comment has been minimized.

Sign in to view

Copy link
@tomidakn

tomidakn Feb 14, 2025

Author Contributor

Oh, I am sorry if this caused a problem. Can you explain more about your trouble?
} // end loop over inner, outer shearing boundaries
} // end shearing box component of ctor
}
Expand Down
10 changes: 8 additions & 2 deletions src/bvals/fc/bvals_shear_emf.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
//========================================================================================

// C headers
#include <float.h>

// C++ headers
#include <algorithm> // min
Expand Down Expand Up @@ -277,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;
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;
}
53 changes: 20 additions & 33 deletions src/bvals/fc/flux_correction_fc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -710,27 +710,27 @@ void FaceCenteredBoundaryVariable::SetFluxBoundarySameLevel(Real *buf,
// store e2 for shearing periodic bcs
for (int k=pmb->ks; k<=pmb->ke+1; k++) {
for (int j=pmb->js; j<=pmb->je; j++)
shear_var_emf_[0].x2e(k,j) =
std::max(shear_var_emf_[0].x2e(k,j), buf[p++]);
shear_var_emf_[0].x2e(k,j)
= std::max(shear_var_emf_[0].x2e(k,j), buf[p++]);
}
// store e3 for shearing periodic bcs
for (int k=pmb->ks; k<=pmb->ke; k++) {
for (int j=pmb->js; j<=pmb->je+1; j++)
shear_var_emf_[0].x3e(k,j) =
std::max(shear_var_emf_[0].x3e(k,j), buf[p++]);
shear_var_emf_[0].x3e(k,j)
= std::max(shear_var_emf_[0].x3e(k,j), buf[p++]);
}
} else if (nb.fid == BoundaryFace::outer_x1) {
// store e2 for shearing periodic bcs
for (int k=pmb->ks; k<=pmb->ke+1; k++) {
for (int j=pmb->js; j<=pmb->je; j++)
shear_var_emf_[1].x2e(k,j) =
std::max(shear_var_emf_[1].x2e(k,j), buf[p++]);
shear_var_emf_[1].x2e(k,j)
= std::max(shear_var_emf_[1].x2e(k,j), buf[p++]);
}
// store e3 for shearing periodic bcs
for (int k=pmb->ks; k<=pmb->ke; k++) {
for (int j=pmb->js; j<=pmb->je+1; j++)
shear_var_emf_[1].x3e(k,j) =
std::max(shear_var_emf_[1].x3e(k,j), buf[p++]);
shear_var_emf_[1].x3e(k,j)
= std::max(shear_var_emf_[1].x3e(k,j), buf[p++]);
}
}
} else {
Expand Down Expand Up @@ -817,10 +817,8 @@ void FaceCenteredBoundaryVariable::SetFluxBoundarySameLevel(Real *buf,
}
} else {
// unpack e2
for (int j=pmb->js; j<=pmb->je; j++) {
e2(k,j,i) = std::max(e2(k,j,i), buf[p++]);
e2(k+1,j,i) = e2(k,j,i);
}
for (int j=pmb->js; j<=pmb->je; j++)
e2(k+1,j,i) = e2(k,j,i) = std::max(e2(k,j,i), buf[p++]);
// unpack e3
for (int j=pmb->js; j<=pmb->je+1; j++)
e3(k,j,i) = std::max(e3(k,j,i), buf[p++]);
Expand All @@ -834,10 +832,8 @@ void FaceCenteredBoundaryVariable::SetFluxBoundarySameLevel(Real *buf,
j = pmb->je + 1;
}
// unpack e1
for (int i=pmb->is; i<=pmb->ie; i++) {
e1(k,j,i) = std::max(e1(k,j,i), buf[p++]);
e1(k+1,j,i) = e1(k,j,i);
}
for (int i=pmb->is; i<=pmb->ie; i++)
e1(k+1,j,i) = e1(k,j,i) = std::max(e1(k,j,i), buf[p++]);
// unpack e3
for (int i=pmb->is; i<=pmb->ie+1; i++)
e3(k,j,i) = std::max(e3(k,j,i), buf[p++]);
Expand All @@ -850,11 +846,9 @@ void FaceCenteredBoundaryVariable::SetFluxBoundarySameLevel(Real *buf,
i = pmb->ie + 1;
}
// unpack e2
e2(k,j,i) = std::max(e2(k,j,i), buf[p++]);
e2(k+1,j,i) = e2(k,j,i);
e2(k+1,j,i) = e2(k,j,i) = std::max(e2(k,j,i), buf[p++]);
// unpack e3
e3(k,j,i) = std::max(e3(k,j,i), buf[p++]);
e3(k,j+1,i) = e3(k,j,i);
e3(k,j+1,i) = e3(k,j,i) = std::max(e3(k,j,i), buf[p++]);
}
} else if (nb.ni.type == NeighborConnect::edge) {
// x1x2 edge (2D and 3D)
Expand Down Expand Up @@ -1058,10 +1052,8 @@ void FaceCenteredBoundaryVariable::SetFluxBoundaryFromFiner(Real *buf,
jl = pmb->js + pmb->block_size.nx2/2;
}
// unpack e2
for (int j=jl; j<=ju; j++) {
e2(k,j,i) = std::max(e2(k,j,i), buf[p++]);
e2(k+1,j,i) = e2(k,j,i);
}
for (int j=jl; j<=ju; j++)
e2(k+1,j,i) = e2(k,j,i) = std::max(e2(k,j,i), buf[p++]);
// unpack e3
for (int j=jl; j<=ju+1; j++)
e3(k,j,i) = std::max(e3(k,j,i), buf[p++]);
Expand All @@ -1079,10 +1071,8 @@ void FaceCenteredBoundaryVariable::SetFluxBoundaryFromFiner(Real *buf,
il = pmb->is + pmb->block_size.nx1/2;
}
// unpack e1
for (int i=il; i<=iu; i++) {
e1(k,j,i) = std::max(e1(k,j,i), buf[p++]);
e1(k+1,j,i) = e1(k,j,i);
}
for (int i=il; i<=iu; i++)
e1(k+1,j,i) = e1(k,j,i) = std::max(e1(k,j,i), buf[p++]);
// unpack e3
for (int i=il; i<=iu+1; i++)
e3(k,j,i) = std::max(e3(k,j,i), buf[p++]);
Expand All @@ -1095,11 +1085,9 @@ void FaceCenteredBoundaryVariable::SetFluxBoundaryFromFiner(Real *buf,
i = pmb->ie + 1;
}
// unpack e2
e2(k,j,i) = std::max(e2(k,j,i), buf[p++]);
e2(k+1,j,i) = e2(k,j,i);
e2(k+1,j,i) = e2(k,j,i) = std::max(e2(k,j,i), buf[p++]);
// unpack e3
e3(k,j,i) = std::max(e3(k,j,i), buf[p++]);
e3(k,j+1,i) = e3(k,j,i);
e3(k,j+1,i) = e3(k,j,i) = std::max(e3(k,j,i), buf[p++]);
}
} else if (nb.ni.type == NeighborConnect::edge) {
if (pmb->block_size.nx3 > 1) { // 3D
Expand Down Expand Up @@ -1547,7 +1535,6 @@ bool FaceCenteredBoundaryVariable::ReceiveFluxCorrection() {
}

if (flag) {
// AverageFluxBoundary();
if (pbval_->num_north_polar_blocks_ > 0)
SetFluxBoundaryFromPolar(flux_north_recv_, pbval_->num_north_polar_blocks_, true);
for (int n = 0; n < pbval_->num_north_polar_blocks_; ++n)
Expand Down

0 comments on commit 1802abe

Please sign in to comment.