From 74365836e29ef0ee6661eb1320a325875d2eaae8 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 3 Jun 2025 13:34:05 +0200 Subject: [PATCH 1/3] fix axis BX2s reconstruction when axis is not used on both sides of the X2 domain (#342) --- src/fluid/boundary/axis.cpp | 53 +++++---------------------------- src/fluid/boundary/axis.hpp | 5 ++-- src/fluid/boundary/boundary.hpp | 46 +++++++++++++++------------- 3 files changed, 37 insertions(+), 67 deletions(-) diff --git a/src/fluid/boundary/axis.cpp b/src/fluid/boundary/axis.cpp index 3ad67f0a0..e447aefc8 100644 --- a/src/fluid/boundary/axis.cpp +++ b/src/fluid/boundary/axis.cpp @@ -358,65 +358,28 @@ void Axis::EnforceAxisBoundary(int side) { // Reconstruct Bx2s taking care of the sides where an axis is lying -void Axis::ReconstructBx2s() { - idfx::pushRegion("Axis::ReconstructBx2s"); +void Axis::RegularizeBX2s() { + idfx::pushRegion("Axis::RegularizeBX2s"); #if DIMENSIONS >= 2 && MHD == YES IdefixArray4D Vs = this->Vs; - IdefixArray3D Ax1=data->A[IDIR]; - IdefixArray3D Ax2=data->A[JDIR]; - IdefixArray3D Ax3=data->A[KDIR]; - int nstart = data->beg[JDIR]-1; - int nend = data->end[JDIR]; - int ntot = data->np_tot[JDIR]; - - bool haveleft = axisLeft; - bool haveright = axisRight; - - if(haveleft) { - // This loop is a copy of ReconstructNormalField, with the proper sign when we cross the axis - idefix_for("Axis::ReconstructBX2sLeft",0,data->np_tot[KDIR],0,data->np_tot[IDIR], - KOKKOS_LAMBDA (int k, int i) { - for(int j = nstart ; j>=0 ; j-- ) { - Vs(BX2s,k,j,i) = 1.0 / Ax2(k,j,i) * ( Ax2(k,j+1,i)*Vs(BX2s,k,j+1,i) - +(D_EXPAND( Ax1(k,j,i+1) * Vs(BX1s,k,j,i+1) - Ax1(k,j,i) * Vs(BX1s,k,j,i) , - , - - ( Ax3(k+1,j,i) * Vs(BX3s,k+1,j,i) - - Ax3(k,j,i) * Vs(BX3s,k,j,i) )))); - } - } - ); - } - if(haveright) { - // This loop is a copy of ReconstructNormalField, with the proper sign when we cross the axis - idefix_for("Axis::ReconstructBX2sRight",0,data->np_tot[KDIR],0,data->np_tot[IDIR], - KOKKOS_LAMBDA (int k, int i) { - for(int j = nend ; jend[JDIR]; int jleft = data->beg[JDIR]; if(isTwoPi) { #ifdef AXIS_BX2S_USE_ATHENA_REGULARISATION - if(haveleft) FixBx2sAxisGhostAverage(left); - if(haveright) FixBx2sAxisGhostAverage(right); + if(axisLeft) FixBx2sAxisGhostAverage(left); + if(axisRight) FixBx2sAxisGhostAverage(right); #else - if(haveleft) FixBx2sAxis(left); - if(haveright) FixBx2sAxis(right); + if(axisLeft) FixBx2sAxis(left); + if(axisRight) FixBx2sAxis(right); #endif } else { + bool haveleft = axisLeft; + bool haveright = axisRight; idefix_for("Axis:BoundaryAvg",0,data->np_tot[KDIR],0,data->np_tot[IDIR], KOKKOS_LAMBDA (int k, int i) { if(haveleft) { diff --git a/src/fluid/boundary/axis.hpp b/src/fluid/boundary/axis.hpp index a1c3054b3..b74697862 100644 --- a/src/fluid/boundary/axis.hpp +++ b/src/fluid/boundary/axis.hpp @@ -29,7 +29,7 @@ class Axis { void RegularizeEMFs(); // Regularize the EMF sitting on the axis void RegularizeCurrent(); // Regularize the currents along the axis void EnforceAxisBoundary(int side); // Enforce the boundary conditions (along X2) - void ReconstructBx2s(); // Reconstruct BX2s in the ghost zone using divB=0 + void RegularizeBX2s(); // Regularize BX2s on the axis void ShowConfig(); @@ -42,6 +42,8 @@ class Axis { void ExchangeMPI(int side); // Function has to be public for GPU, but its technically // a private function + bool haveLeftAxis() const { return axisLeft; } ///< Check if the axis is on the left + bool haveRightAxis() const { return axisRight; } ///< Check if the axis is on the right private: bool isTwoPi = false; @@ -147,7 +149,6 @@ Axis::Axis(Boundary *boundary) { Kokkos::deep_copy(symmetryVc, symmetryVcHost); if constexpr(Phys::mhd) { - idfx::cout << "Phys MHD" << std::endl; symmetryVs = IdefixArray1D("Axis:SymmetryVs",DIMENSIONS); IdefixArray1D::HostMirror symmetryVsHost = Kokkos::create_mirror_view(symmetryVs); // Init the array diff --git a/src/fluid/boundary/boundary.hpp b/src/fluid/boundary/boundary.hpp index 86909936f..108793e1d 100644 --- a/src/fluid/boundary/boundary.hpp +++ b/src/fluid/boundary/boundary.hpp @@ -112,6 +112,8 @@ class Boundary { IdefixArray4D Vs; ///< reference to face-centered array that we should sync std::unique_ptr axis; ///< Axis object, initialised if needed. bool haveAxis{false}; + bool haveLeftAxis{false}; ///< True if the left boundary is an axis + bool haveRightAxis{false}; ///< True if the right boundary is an axis private: friend class Axis; @@ -137,6 +139,8 @@ Boundary::Boundary(Fluid* fluid) { if(data->haveAxis) { this->axis = std::make_unique(this); this->haveAxis = true; + this->haveLeftAxis = axis->haveLeftAxis(); + this->haveRightAxis = axis->haveRightAxis(); } @@ -464,30 +468,32 @@ void Boundary::ReconstructNormalField(int dir) { if(dir==JDIR) { nstart = data->beg[JDIR]-1; nend = data->end[JDIR]; - if(!this->haveAxis) { - idefix_for("ReconstructBX2s",0,data->np_tot[KDIR],0,data->np_tot[IDIR], - KOKKOS_LAMBDA (int k, int i) { - if(reconstructLeft) { - for(int j = nstart ; j>=0 ; j-- ) { - Vs(BX2s,k,j,i) = 1.0 / Ax2(k,j,i) * ( Ax2(k,j+1,i)*Vs(BX2s,k,j+1,i) - +(D_EXPAND( Ax1(k,j,i+1) * Vs(BX1s,k,j,i+1) - Ax1(k,j,i) * Vs(BX1s,k,j,i) , - , - + Ax3(k+1,j,i) * Vs(BX3s,k+1,j,i) - Ax3(k,j,i) * Vs(BX3s,k,j,i) ))); - } + const int signLeft = haveLeftAxis ? -1 : 1; // left axis is in the negative direction + const int signRight = haveRightAxis ? -1 : 1; // right axis is in the negative direction + + idefix_for("ReconstructBX2s",0,data->np_tot[KDIR],0,data->np_tot[IDIR], + KOKKOS_LAMBDA (int k, int i) { + if(reconstructLeft) { + for(int j = nstart ; j>=0 ; j-- ) { + Vs(BX2s,k,j,i) = 1.0 / Ax2(k,j,i) * ( Ax2(k,j+1,i)*Vs(BX2s,k,j+1,i) + +(D_EXPAND( Ax1(k,j,i+1) * Vs(BX1s,k,j,i+1) - Ax1(k,j,i) * Vs(BX1s,k,j,i) , + , + + signLeft*(Ax3(k+1,j,i) * Vs(BX3s,k+1,j,i) - Ax3(k,j,i) * Vs(BX3s,k,j,i) )))); } - if(reconstructRight) { - for(int j = nend ; jReconstructBx2s(); + axis->RegularizeBX2s(); } } #endif From 0c16b8cdf7242edf8b4f50df304e54a05013b0dc Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 3 Jun 2025 14:01:03 +0200 Subject: [PATCH 2/3] fix reflective with MHD and DIMENSIONS < 3 (#343) --- src/fluid/boundary/boundary.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fluid/boundary/boundary.hpp b/src/fluid/boundary/boundary.hpp index 108793e1d..8be40a3fb 100644 --- a/src/fluid/boundary/boundary.hpp +++ b/src/fluid/boundary/boundary.hpp @@ -692,7 +692,7 @@ void Boundary::EnforceReflective(int dir, BoundarySide side ) { const int jref = (dir==JDIR) ? 2*(jghost + side*nxj) - j - 1 : j; const int kref = (dir==KDIR) ? 2*(kghost + side*nxk) - k - 1 : k; - const int sign = (n == VX1+dir) ? -1.0 : 1.0; + const int sign = (n == VX1+dir || (n >= BX1 && n != BX1+dir)) ? -1.0 : 1.0; Vc(n,k,j,i) = sign * Vc(n,kref,jref,iref); }); From b35b555985cc66db7981df1847a96527d428feac Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 5 Jun 2025 10:27:57 +0200 Subject: [PATCH 3/3] fix warning on cuda --- src/fluid/boundary/boundary.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/fluid/boundary/boundary.hpp b/src/fluid/boundary/boundary.hpp index 8be40a3fb..897c1ce10 100644 --- a/src/fluid/boundary/boundary.hpp +++ b/src/fluid/boundary/boundary.hpp @@ -468,8 +468,10 @@ void Boundary::ReconstructNormalField(int dir) { if(dir==JDIR) { nstart = data->beg[JDIR]-1; nend = data->end[JDIR]; - const int signLeft = haveLeftAxis ? -1 : 1; // left axis is in the negative direction - const int signRight = haveRightAxis ? -1 : 1; // right axis is in the negative direction + #if DIMENSIONS == 3 + const int signLeft = haveLeftAxis ? -1 : 1; // left axis is in the negative direction + const int signRight = haveRightAxis ? -1 : 1; // right axis is in the negative direction + #endif idefix_for("ReconstructBX2s",0,data->np_tot[KDIR],0,data->np_tot[IDIR], KOKKOS_LAMBDA (int k, int i) {