From 6bd24251b8699bd505e5323c86f90cd77f02d964 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 23 Jun 2026 21:54:54 +0200 Subject: [PATCH 01/11] add mean field documentation --- doc/source/input_file.rst | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/doc/source/input_file.rst b/doc/source/input_file.rst index e07edc1..d34846c 100644 --- a/doc/source/input_file.rst +++ b/doc/source/input_file.rst @@ -159,7 +159,11 @@ The ``InitFlow`` section defines the initial conditions of the flow. Several ent | shear_layer | float, float | | (optional) Create a shear layer: :math:`v_x=-v_0` if :math:`y Date: Tue, 23 Jun 2026 21:55:11 +0200 Subject: [PATCH 02/11] clean up hydro class --- src/rightHandSide/hydro.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/rightHandSide/hydro.hpp b/src/rightHandSide/hydro.hpp index a51f8ef..8ba94b8 100644 --- a/src/rightHandSide/hydro.hpp +++ b/src/rightHandSide/hydro.hpp @@ -305,7 +305,6 @@ template void Hydro::PostStage(Field>& fldin, real t) { if constexpr(Shear::isEnabled) { if(this->shear.NeedRemap(t)) { - astra::cout << "Hydro rhs: remapping fields at t=" << t << std::endl; this->shear.Remap(t, fldin["vx1"]); this->shear.Remap(t, fldin["vx2"]); this->shear.Remap(t, fldin["vx3"]); From cb93ee1568344ba6fec236876273064b230a5848 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 23 Jun 2026 21:55:27 +0200 Subject: [PATCH 03/11] duplicated source term in mhd --- src/rightHandSide/mhd.hpp | 26 -------------------------- 1 file changed, 26 deletions(-) diff --git a/src/rightHandSide/mhd.hpp b/src/rightHandSide/mhd.hpp index f3badde..84deab6 100644 --- a/src/rightHandSide/mhd.hpp +++ b/src/rightHandSide/mhd.hpp @@ -188,32 +188,6 @@ void Mhd::ExplicitStep(Field>& fldin, FieldOmega; - real S = this->shear.shearRate; - astra_for("hydro_source_terms", 0,npf[IDIR],0,npf[JDIR],0,npf[KDIR], - KOKKOS_LAMBDA(int64_t i, int64_t j, int64_t k) { - dvx1(i,j,k) += 2.0*Omega*vx2(i,j,k); - dvx2(i,j,k) += -(2.0*Omega - S)*vx1(i,j,k); - }); - } - // Pressure term - astra_for("hydro_pressure", 0,npf[IDIR],0,npf[JDIR],0,npf[KDIR], - KOKKOS_LAMBDA(int64_t i, int64_t j, int64_t k) { - const real kx1t = shear.kx1t(kx1(i),kx2(j),kx3(k)); - const real kx2t = shear.kx2t(kx1(i),kx2(j),kx3(k)); - const real kx3t = shear.kx3t(kx1(i),kx2(j),kx3(k)); - const real k2t = kx1t*kx1t + kx2t*kx2t + kx3t*kx3t; - if(k2t > 0.0) { - complex kv_dot_v = kx1t*dvx1(i,j,k) + kx2t*dvx2(i,j,k) + kx3t*dvx3(i,j,k); - kv_dot_v += shear.shearRate * kx2(j) * vx1(i,j,k); // Shear contribution = dk/dt.v - dvx1(i,j,k) -= kv_dot_v*kx1t/k2t; - dvx2(i,j,k) -= kv_dot_v*kx2t/k2t; - dvx3(i,j,k) -= kv_dot_v*kx3t/k2t; - } - }); - // Magnetic field induction // vxB in real space astra_for("mhd_emf", 0,npr[IDIR],0,npr[JDIR],0,npr[KDIR], From 8af64fc1cc08be0e48173d2825268dbbbc5991f8 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 23 Jun 2026 21:55:47 +0200 Subject: [PATCH 04/11] apply 2/3 antialiasing to remap procedure --- src/shear/linearshear.hpp | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/shear/linearshear.hpp b/src/shear/linearshear.hpp index a8e5e99..edef38d 100644 --- a/src/shear/linearshear.hpp +++ b/src/shear/linearshear.hpp @@ -58,6 +58,13 @@ class LinearShear : public NoShear { int nx_glob = this->grid->npf_glob[IDIR]; int ny_glob = this->grid->npf_glob[JDIR]; + + real kx1max = this->grid->kmax[IDIR]; + real kx2max = this->grid->kmax[JDIR]; + real kx3max = this->grid->kmax[KDIR]; + auto kx1 = this->grid->kx[IDIR]; + auto kx2 = this->grid->kx[JDIR]; + auto kx3 = this->grid->kx[KDIR]; #ifdef WITH_MPI Array3D transposed("transposed", this->grid->npf_glob[JDIR]/astra::psize, this->grid->npf_glob[IDIR], this->grid->npf_glob[KDIR]); Transpose transpose(this->grid->npf); @@ -76,7 +83,10 @@ class LinearShear : public NoShear { // Check if mode goes out of bounds if(nxtarget > -nx_glob/2 && nxtarget <= nx_glob/2) { const int inew = (nxtarget + nx_glob) % nx_glob; - temp(j,inew,k) = transposed(j,i,k); + complex mask = (std::fabs(kx1(inew))< 2./3*kx1max + && std::fabs(kx2(j))< 2./3*kx2max + && std::fabs(kx3(k))< 2./3*kx3max) ? 1.0 : 0.0; + temp(j,inew,k) = mask*transposed(j,i,k); } } ); @@ -95,7 +105,10 @@ class LinearShear : public NoShear { // Check if mode goes out of bounds if(nxtarget > -nx_glob/2 && nxtarget <= nx_glob/2) { const int inew = (nxtarget + nx_glob) % nx_glob; - temp(inew,j,k) = field(i,j,k); + complex mask = (std::fabs(kx1(inew))< 2./3*kx1max + && std::fabs(kx2(j))< 2./3*kx2max + && std::fabs(kx3(k))< 2./3*kx3max) ? 1.0 : 0.0; + temp(inew,j,k) = mask*field(i,j,k); } } ); From 10e0ba0a762aae0f2e85a1ffc41e0c28892991fb Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 23 Jun 2026 21:56:27 +0200 Subject: [PATCH 05/11] change cfl and resolution for ci hydro test --- test/shearing_wave/astra.ini | 4 ++-- test/shearing_wave/validate.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/shearing_wave/astra.ini b/test/shearing_wave/astra.ini index 0206b33..71c1426 100644 --- a/test/shearing_wave/astra.ini +++ b/test/shearing_wave/astra.ini @@ -1,5 +1,5 @@ [Grid] -X1-grid -0.5 32 0.5 +X1-grid -0.5 64 0.5 X2-grid -0.5 32 0.5 X3-grid -0.5 32 0.5 @@ -8,7 +8,7 @@ script shearing_wave_test [TimeIntegrator] method rk3 -cfl 0.1 +cfl 0.9 tstop 10.0 [Physics] diff --git a/test/shearing_wave/validate.py b/test/shearing_wave/validate.py index 81a249b..e7352d9 100755 --- a/test/shearing_wave/validate.py +++ b/test/shearing_wave/validate.py @@ -107,7 +107,7 @@ def rhs(t, y, Omega, q, k0x, k0y, k0z): err=np.mean(error) print("Error=",err) -if(err<8e-4): +if(err<7e-3): print("SUCCESS") sys.exit(0) else: From 188aa5ad398c5bc2c6afcd3de96702540c29cb91 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 23 Jun 2026 21:56:59 +0200 Subject: [PATCH 06/11] change resolution for compressible hydro test --- test/shearing_wave_compressible/astra.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/shearing_wave_compressible/astra.ini b/test/shearing_wave_compressible/astra.ini index 6caed90..d9494d6 100644 --- a/test/shearing_wave_compressible/astra.ini +++ b/test/shearing_wave_compressible/astra.ini @@ -1,5 +1,5 @@ [Grid] -X1-grid -0.5 32 0.5 +X1-grid -0.5 64 0.5 X2-grid -0.5 32 0.5 X3-grid -0.5 32 0.5 From fb4beff8a30dad3231934e326de19d28039d9e0e Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 23 Jun 2026 21:58:22 +0200 Subject: [PATCH 07/11] add mhd shearing wave test --- test/shearing_wave_mhd/astra.ini | 28 ++++ test/shearing_wave_mhd/shearing_wave_test.py | 28 ++++ test/shearing_wave_mhd/validate.py | 143 +++++++++++++++++++ 3 files changed, 199 insertions(+) create mode 100644 test/shearing_wave_mhd/astra.ini create mode 100644 test/shearing_wave_mhd/shearing_wave_test.py create mode 100755 test/shearing_wave_mhd/validate.py diff --git a/test/shearing_wave_mhd/astra.ini b/test/shearing_wave_mhd/astra.ini new file mode 100644 index 0000000..7432a9b --- /dev/null +++ b/test/shearing_wave_mhd/astra.ini @@ -0,0 +1,28 @@ +[Grid] +X1-grid -0.5 64 0.5 +X2-grid -0.5 32 0.5 +X3-grid -0.5 32 0.5 + +[Python] +script shearing_wave_test + +[TimeIntegrator] +method rk3 +cfl 0.9 +tstop 10.0 + +[Physics] +rhs mhd +resistivity 0 +viscosity 0 +omega 1.0 +shear_type linear +shear_rate 1.5 + +[InitFlow] +python init_flow fourier +mean_field 0.0 0.02 0.05 # Mean magnetic field in each direction + +[Output] +log 10 +python 0.1 output fourier diff --git a/test/shearing_wave_mhd/shearing_wave_test.py b/test/shearing_wave_mhd/shearing_wave_test.py new file mode 100644 index 0000000..c802dc2 --- /dev/null +++ b/test/shearing_wave_mhd/shearing_wave_test.py @@ -0,0 +1,28 @@ +from astrapy import * +import numpy as np +#import matplotlib.pyplot as plt + +# The output function +# the only argument is dataBlockHost python object, wrapping a dataBlockHost Idefix object +def output(grid, field, t, n): + # only process #0 performs the output + shear = 1.5 + if prank==0: + mode = "a" + if n==0: + mode = "w" + f = open("mode_amplitude.dat",mode) + kx = int(t*shear+0.5) + if n==0: + f.write("t vx vy vz bx by bz\n") + + f.write("%e %e %e %e %e %e %e\n" % (t, np.real(field["vx1"][kx,1,4]), np.real(field["vx2"][kx,1,4]), np.real(field["vx3"][kx,1,4]), -np.imag(field["bx1"][kx,1,4]), -np.imag(field["bx2"][kx,1,4]), -np.imag(field["bx3"][kx,1,4]))) + f.close() + + +def init_flow(grid, field): + if prank==0: + field["vx1"][0,1,4] = 1.0e-2*grid.npr_glob[IDIR]*grid.npr_glob[JDIR]*grid.npr_glob[KDIR] + + # Field amplitude + #field["vx1"] = 0*x diff --git a/test/shearing_wave_mhd/validate.py b/test/shearing_wave_mhd/validate.py new file mode 100755 index 0000000..c7fd0f1 --- /dev/null +++ b/test/shearing_wave_mhd/validate.py @@ -0,0 +1,143 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Jun 21 15:42:19 2021 + +@author: lesurg +""" +import sys +import numpy as np +import matplotlib.pyplot as plt +from scipy.integrate import solve_ivp +import argparse + +#compute theoretical solution (from Balbus & Hawley 2006, using notations from Lesur 2021) +def rhs(t, y, Omega, q, B0y, B0z, k0x, k0y, k0z): + kx = k0x + q*Omega * t * k0y + ky = k0y + kz = k0z + + k2 = kx*kx+ky*ky+kz*kz + vx = y[0] + vy = y[1] + vz = y[2] + bx = y[3] + by = y[4] + bz = y[5] + + gxx = kx*kx/k2 + gxy = kx*ky/k2 + gyy = ky*ky/k2 + gyz = ky*kz/k2 + gxz = kx*kz/k2 + #gzz = kz*kz/k2 + + kdotB = ky*B0y+kz*B0z + + #we assume b=1j*db hence db=-1j*b + + dvx = kdotB*bx + 2*Omega*vy*(1-gxx) + 2*(1-q)*Omega*gxy*vx + dvy = kdotB*by - q*Omega*vx*gyy - (2-q)*Omega*(1-gyy)*vx - 2*Omega*vy*gxy + dvz= kdotB*bz +2*(1-q)*Omega*vx*gyz-2*Omega*vy*gxz + + dbx = - kdotB*vx + dby = - kdotB*vy - q*Omega*bx + dbz = - kdotB*vz + + return np.asarray([dvx,dvy,dvz,dbx,dby,dbz]) + +parser = argparse.ArgumentParser() +parser.add_argument("-noplot", + default=False, + help="disable plotting", + action="store_true") + + +args, unknown=parser.parse_known_args() + + +#initial condition: vr=1, rest is 0, mode initial is nx=0, ny=1, nz=4) +y=solve_ivp(rhs, [0,30], [1, 0, 0, 0, 0, 0], args=(1, 1.5, 0.02, 0.05, 0, 2.0*np.pi, 8*np.pi), dense_output=True) + + + +# read timevol file +rep = "./" +filename = "mode_amplitude.dat" + +fid=open(rep+filename,"r") + +gamma=5/3 +# read the first line to get data names +varnames=fid.readline().split() +fid.close() + +# load the bulk of the file +data=np.loadtxt(rep+filename,skiprows=1) + +# store this in our data structure +V={} + +i=0 +for name in varnames: + V[name]=data[:,i] + i=i+1 + +# velocity normalisation +v0=V['vx'][0] +# compute L2 error norm +time = V["t"] +error = (V['vx']/v0-y.sol(time)[0,:])**2 +error = error + (V['vy']/v0-y.sol(time)[1,:])**2 +error = error + (V['vz']/v0-y.sol(time)[2,:])**2 +error = error + (V['bx']/v0-y.sol(time)[3,:])**2 +error = error + (V['by']/v0-y.sol(time)[4,:])**2 +error = error + (V['bz']/v0-y.sol(time)[5,:])**2 +error = np.sqrt(error) + + + +if(not args.noplot): + # Comparison with exact solution + plt.rc('text', usetex=True) + plt.rc('font', family='serif') + plt.rc('font', size=16) + plt.close('all') + plt.figure(1) + plt.plot(V["t"],V['vx']/v0,'r-',label=r'$u_{R}$') + plt.plot(V["t"],y.sol(V["t"])[0,:],'r--') + plt.plot(V["t"],V['vy']/v0,'b-',label=r'$u_{\varphi}$') + plt.plot(V["t"],y.sol(V["t"])[1,:],'b--') + plt.plot(V["t"],V['vz']/v0,'g-',label=r'$u_{z}$') + plt.plot(V["t"],y.sol(V["t"])[2,:],'g--') + 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) + plt.xlim([0,10]) + #plt.ylim([1e-3,1]) + + plt.ioff() + plt.show() + +err=np.mean(error) +print("Error=",err) + +if(err<8e-3): + print("SUCCESS") + sys.exit(0) +else: + print("FAILED") + sys.exit(1) From c95bf8e28410aed4a5e13b3ce48fcb2054e7a145 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 23 Jun 2026 21:59:33 +0200 Subject: [PATCH 08/11] add compressible hydro and mhd tests to CI --- .github/workflows/astra-ci-jobs.yml | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/.github/workflows/astra-ci-jobs.yml b/.github/workflows/astra-ci-jobs.yml index 02b87aa..9f3dfd8 100644 --- a/.github/workflows/astra-ci-jobs.yml +++ b/.github/workflows/astra-ci-jobs.yml @@ -46,6 +46,16 @@ jobs: source .venv/bin/activate cd job-serial python3 ../test/test.py run shearing_wave $TEST_OPTIONS + - name: MHD tests + run: | + source .venv/bin/activate + cd job-serial + python3 ../test/test.py run shearing_wave_mhd $TEST_OPTIONS + - name: Compressible hydro tests + run: | + source .venv/bin/activate + cd job-serial + python3 ../test/test.py run shearing_wave_compressible $TEST_OPTIONS Parallel: runs-on: self-hosted @@ -79,3 +89,13 @@ jobs: source .venv/bin/activate cd job-parallel python3 ../test/test.py run shearing_wave $TEST_OPTIONS -mpi + - name: MHD tests + run: | + source .venv/bin/activate + cd job-parallel + python3 ../test/test.py run shearing_wave_mhd $TEST_OPTIONS -mpi + - name: Compressible hydro tests + run: | + source .venv/bin/activate + cd job-parallel + python3 ../test/test.py run shearing_wave_compressible $TEST_OPTIONS -mpi From 45e13e5f57dcc081e2dcbcb9d0a893c6491986b0 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 23 Jun 2026 22:09:37 +0200 Subject: [PATCH 09/11] don't make vtks for compressible shearing wave --- test/shearing_wave_compressible/astra.ini | 1 - 1 file changed, 1 deletion(-) diff --git a/test/shearing_wave_compressible/astra.ini b/test/shearing_wave_compressible/astra.ini index d9494d6..8b3be6a 100644 --- a/test/shearing_wave_compressible/astra.ini +++ b/test/shearing_wave_compressible/astra.ini @@ -26,5 +26,4 @@ python init_flow fourier [Output] log 10 -vtk 0.1 python 0.1 output fourier From 37c6de69710fe1f14108b4905dd44ede9afe5a2a Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Tue, 23 Jun 2026 22:55:16 +0200 Subject: [PATCH 10/11] apply mask after transposition when remapping --- src/shear/linearshear.hpp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/shear/linearshear.hpp b/src/shear/linearshear.hpp index edef38d..330583b 100644 --- a/src/shear/linearshear.hpp +++ b/src/shear/linearshear.hpp @@ -83,9 +83,6 @@ class LinearShear : public NoShear { // Check if mode goes out of bounds if(nxtarget > -nx_glob/2 && nxtarget <= nx_glob/2) { const int inew = (nxtarget + nx_glob) % nx_glob; - complex mask = (std::fabs(kx1(inew))< 2./3*kx1max - && std::fabs(kx2(j))< 2./3*kx2max - && std::fabs(kx3(k))< 2./3*kx3max) ? 1.0 : 0.0; temp(j,inew,k) = mask*transposed(j,i,k); } } @@ -105,15 +102,22 @@ class LinearShear : public NoShear { // Check if mode goes out of bounds if(nxtarget > -nx_glob/2 && nxtarget <= nx_glob/2) { const int inew = (nxtarget + nx_glob) % nx_glob; - complex mask = (std::fabs(kx1(inew))< 2./3*kx1max - && std::fabs(kx2(j))< 2./3*kx2max - && std::fabs(kx3(k))< 2./3*kx3max) ? 1.0 : 0.0; - temp(inew,j,k) = mask*field(i,j,k); + temp(inew,j,k) = field(i,j,k); } } ); Kokkos::deep_copy(field, temp); #endif + + // Apply 2/3 de-aliasing rule + astra_for("mask_remap", 0,field.extent(0),0,field.extent(1),0,field.extent(2), + KOKKOS_LAMBDA (const int i, const int j, const int k) { + complex mask = (std::fabs(kx1(i))< 2./3*kx1max + && std::fabs(kx2(j))< 2./3*kx2max + && std::fabs(kx3(k))< 2./3*kx3max) ? 1.0 : 0.0; + field(i,j,k) *= mask; + } + ); // Done astra::popRegion(); } From 3a7bb5c822e9e27f95ac99437b918c8d0ee73fef Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 24 Jun 2026 07:20:22 +0200 Subject: [PATCH 11/11] clean up linearshear --- src/shear/linearshear.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/shear/linearshear.hpp b/src/shear/linearshear.hpp index 330583b..948ba4b 100644 --- a/src/shear/linearshear.hpp +++ b/src/shear/linearshear.hpp @@ -83,7 +83,7 @@ class LinearShear : public NoShear { // Check if mode goes out of bounds if(nxtarget > -nx_glob/2 && nxtarget <= nx_glob/2) { const int inew = (nxtarget + nx_glob) % nx_glob; - temp(j,inew,k) = mask*transposed(j,i,k); + temp(j,inew,k) = transposed(j,i,k); } } );