Mhd shearing box

.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

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<y_0`, else :math:`v_x=v_0`               |
 |                        |                    | | The first parameter is :math:`y_0` and the second parameter is :math:`v_0`.                             |
 +------------------------+--------------------+-----------------------------------------------------------------------------------------------------------+
-| mean_flow              | float, float, float|  (optional) Create a mean flow. The three parameters are the three components of the mean velocity        |
+| mean_flow              | float, float, float|  (optional) Adds a mean flow. The three parameters are the three components of the mean velocity          |
++------------------------+--------------------+-----------------------------------------------------------------------------------------------------------+
+| mean_field             | float, float, float| | (optional) Adds a constant magnetic field. The three parameters are the three components of field       |
+|                        |                    | | :math:`\mathbf{B}=(B_x, B_y, B_z)`.                                                                     |
+|                        |                    | | NB: this entry is only available when the rhs is ``mhd``.                                               |
 +------------------------+--------------------+-----------------------------------------------------------------------------------------------------------+
 | python                 | string, (string)   | | 1st parameter: Name of the Python function to call to initialize the flow in the script provided in the |
 |                        |                    | | [python] block of the input file.                   ,                                                   |

src/rightHandSide/hydro.hpp

@@ -305,7 +305,6 @@ template <typename Shear>
 void Hydro<Shear>::PostStage(Field<Array3D<complex>>& 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"]);

src/rightHandSide/mhd.hpp

@@ -188,32 +188,6 @@ void Mhd<Shear>::ExplicitStep(Field<Array3D<complex>>& fldin, Field<Array3D<comp
       dvx3(i,j,k) -= Kokkos::complex(0.0,1.0)*(kx1t*wf13(i,j,k)+kx2t*wf23(i,j,k)+kx3t*wf33(i,j,k))*mask;
   });
 
-  // Source terms
-  if(haveSourceTerm) {
-    real Omega = this->Omega;
-    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],

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<complex> transposed("transposed", this->grid->npf_glob[JDIR]/astra::psize, this->grid->npf_glob[IDIR], this->grid->npf_glob[KDIR]);
       Transpose<complex> transpose(this->grid->npf);
@@ -101,6 +108,16 @@ class LinearShear : public NoShear {
       );
       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();
   }

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]

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:

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
 
@@ -26,5 +26,4 @@ python    init_flow  fourier
 
 [Output]
 log       10
-vtk       0.1
 python    0.1  output  fourier

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

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

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)