From 45d5e9fb1652680884488faa8b8ffb79f54589ca Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Mon, 14 Sep 2026 18:14:48 +0200 Subject: [PATCH 1/2] ensure shearing box test also analyse B --- test/MHD/ShearingBox/analysis.cpp | 11 ++++++----- test/MHD/ShearingBox/analysis.hpp | 3 ++- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/test/MHD/ShearingBox/analysis.cpp b/test/MHD/ShearingBox/analysis.cpp index 95102b993..7b93d54f7 100644 --- a/test/MHD/ShearingBox/analysis.cpp +++ b/test/MHD/ShearingBox/analysis.cpp @@ -58,7 +58,8 @@ double Analysis::ShwaveAmplitude(const int field, const int nx, const int ny, const int nz, - const real t) + const real t, + const real phaseShift = 0.0) /* * compute the weighted average: int dphi dz rho *infield/int dphi dz rho * @@ -74,7 +75,7 @@ double Analysis::ShwaveAmplitude(const int field, real z = d->x[KDIR](k); real wave = sin(2.0*M_PI*( (nx-ny*shear*t)*x + ny*y - + nz*z)); + + nz*z + 2*phaseShift*M_PI)); q += wave*d->Vc(field,k,j,i); } } @@ -159,9 +160,9 @@ void Analysis::PerformAnalysis(DataBlock &data) { WriteField(ShwaveAmplitude(VX1, 0, 1, 2, data.t) ); WriteField(ShwaveAmplitude(VX2, 0, 1, 2, data.t) ); WriteField(ShwaveAmplitude(VX3, 0, 1, 2, data.t) ); - WriteField(ShwaveAmplitude(BX1, 0, 1, 2, data.t) ); - WriteField(ShwaveAmplitude(BX2, 0, 1, 2, data.t) ); - WriteField(ShwaveAmplitude(BX3, 0, 1, 2, data.t) ); + WriteField(ShwaveAmplitude(BX1, 0, 1, 2, data.t, -0.5)); + WriteField(ShwaveAmplitude(BX2, 0, 1, 2, data.t, -0.5)); + WriteField(ShwaveAmplitude(BX3, 0, 1, 2, data.t, -0.5)); if(idfx::prank==0) { diff --git a/test/MHD/ShearingBox/analysis.hpp b/test/MHD/ShearingBox/analysis.hpp index 8b65c4317..3a0633f37 100644 --- a/test/MHD/ShearingBox/analysis.hpp +++ b/test/MHD/ShearingBox/analysis.hpp @@ -25,7 +25,8 @@ class Analysis { const int, const int, const int, - const real); + const real, + const real); DataBlockHost *d; Grid *grid; From 713aaaf8afe7a1247e8da6eb87a51fd28ccd2e37 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 17 Sep 2026 09:13:41 +0200 Subject: [PATCH 2/2] add B field in the validation of the shearing box test. Adjust tolerance accordingly. --- test/MHD/ShearingBox/analysis.cpp | 8 ++++---- test/MHD/ShearingBox/python/testidefix.py | 17 +++++++++++++++-- 2 files changed, 19 insertions(+), 6 deletions(-) diff --git a/test/MHD/ShearingBox/analysis.cpp b/test/MHD/ShearingBox/analysis.cpp index 7b93d54f7..c227c3fa8 100644 --- a/test/MHD/ShearingBox/analysis.cpp +++ b/test/MHD/ShearingBox/analysis.cpp @@ -75,7 +75,7 @@ double Analysis::ShwaveAmplitude(const int field, real z = d->x[KDIR](k); real wave = sin(2.0*M_PI*( (nx-ny*shear*t)*x + ny*y - + nz*z + 2*phaseShift*M_PI)); + + nz*z ) + 2*phaseShift*M_PI); q += wave*d->Vc(field,k,j,i); } } @@ -160,9 +160,9 @@ void Analysis::PerformAnalysis(DataBlock &data) { WriteField(ShwaveAmplitude(VX1, 0, 1, 2, data.t) ); WriteField(ShwaveAmplitude(VX2, 0, 1, 2, data.t) ); WriteField(ShwaveAmplitude(VX3, 0, 1, 2, data.t) ); - WriteField(ShwaveAmplitude(BX1, 0, 1, 2, data.t, -0.5)); - WriteField(ShwaveAmplitude(BX2, 0, 1, 2, data.t, -0.5)); - WriteField(ShwaveAmplitude(BX3, 0, 1, 2, data.t, -0.5)); + WriteField(ShwaveAmplitude(BX1, 0, 1, 2, data.t, -0.25)); + WriteField(ShwaveAmplitude(BX2, 0, 1, 2, data.t, -0.25)); + WriteField(ShwaveAmplitude(BX3, 0, 1, 2, data.t, -0.25)); if(idfx::prank==0) { diff --git a/test/MHD/ShearingBox/python/testidefix.py b/test/MHD/ShearingBox/python/testidefix.py index 6f75c6af7..559c6c6ff 100755 --- a/test/MHD/ShearingBox/python/testidefix.py +++ b/test/MHD/ShearingBox/python/testidefix.py @@ -64,7 +64,7 @@ def rhs(t, y, Omega, q, B0y, B0z, k0x, k0y, k0z): args, unknown = parser.parse_known_args() -# initial condition: vr=1, rest is 0, mode initial is nx=0, ny=1, nz=5) +# initial condition: vr=1, rest is 0, mode initial is nx=0, ny=1, nz=2) y = solve_ivp( rhs, [0, 30], @@ -103,6 +103,9 @@ def rhs(t, y, Omega, q, B0y, B0z, k0x, k0y, k0z): (V["vx"] / v0 - y.sol(V["t"])[0, :]) ** 2 + (V["vy"] / v0 - y.sol(V["t"])[1, :]) ** 2 + (V["vz"] / v0 - y.sol(V["t"])[2, :]) ** 2 + + (V["bx"] / v0 - y.sol(V["t"])[3, :]) ** 2 + + (V["by"] / v0 - y.sol(V["t"])[4, :]) ** 2 + + (V["bz"] / v0 - y.sol(V["t"])[5, :]) ** 2 ) @@ -122,6 +125,16 @@ def rhs(t, y, Omega, q, B0y, B0z, k0x, k0y, k0z): plt.legend() plt.xlabel("t") + plt.figure(2) + plt.plot(V["t"], V["bx"] / v0, "r-", label=r"$b_{R}$") + plt.plot(V["t"], y.sol(V["t"])[3, :], "r--") + plt.plot(V["t"], V["by"] / v0, "b-", label=r"$b_{\varphi}$") + plt.plot(V["t"], y.sol(V["t"])[4, :], "b--") + plt.plot(V["t"], V["bz"] / v0, "g-", label=r"$b_{z}$") + plt.plot(V["t"], y.sol(V["t"])[5, :], "g--") + plt.legend() + plt.xlabel("t") + # plot error plt.figure() plt.semilogy(V["t"], error) @@ -134,7 +147,7 @@ def rhs(t, y, Omega, q, B0y, B0z, k0x, k0y, k0z): err = np.mean(error) print("Error=", err) -if err < 0.03: +if err < 0.04: print("SUCCESS") sys.exit(0) else: