Fix parabolic MHD+coarsening (#233)

* adding a new test for grid Coarsening

* fix currents used with grid Coarsening

.github/workflows/idefix-ci.yml

@@ -123,6 +123,10 @@ jobs:
         run: |
           cd $IDEFIX_DIR/test/MHD/ResistiveAlfvenWave
           ./testme.py -all $TESTME_OPTIONS
+      - name: Grid coarsening diffusion
+        run: |
+          cd $IDEFIX_DIR/test/MHD/Coarsening
+          ./testme.py -all $TESTME_OPTIONS
       - name: Hall whistler waves
         run: |
           cd $IDEFIX_DIR/test/MHD/HallWhistler

src/fluid/calcCurrent.hpp

@@ -31,6 +31,14 @@ void Fluid<Phys>::CalcCurrent() {
   IdefixArray1D<real> sinx2m = data->sinx2m;
   #endif
 
+  bool haveGridCoarsening = data->haveGridCoarsening != GridCoarsening::disabled;
+  bool haveGridCoarseningX1 = data->coarseningDirection[IDIR];
+  bool haveGridCoarseningX2 = data->coarseningDirection[JDIR];
+  bool haveGridCoarseningX3 = data->coarseningDirection[KDIR];
+  auto coarseningLevelX1 = data->coarseningLevel[IDIR];
+  auto coarseningLevelX2 = data->coarseningLevel[JDIR];
+  auto coarseningLevelX3 = data->coarseningLevel[KDIR];
+
   idefix_for("CalcCurrent",
              KOFFSET,data->np_tot[KDIR],
              JOFFSET,data->np_tot[JDIR],
@@ -127,7 +135,28 @@ void Fluid<Phys>::CalcCurrent() {
 
 #endif
 
-      // Compute actual current
+    // Correct spacing for grid coarsening
+    if(haveGridCoarsening) {
+      if(haveGridCoarseningX1) {
+        const int factor =  1 << (coarseningLevelX1(k,j) - 1);
+        d12 /= factor;
+        d13 /= factor;
+      }
+
+      if(haveGridCoarseningX2) {
+        const int factor =  1 << (coarseningLevelX2(k,i) - 1);
+        d21 /= factor;
+        d23 /= factor;
+      }
+
+      if(haveGridCoarseningX3) {
+        const int factor =  1 << (coarseningLevelX3(j,i) - 1);
+        d31 /= factor;
+        d32 /= factor;
+      }
+    }
+
+    // Compute actual current
 
 #if COMPONENTS == 3
   #if DIMENSIONS == 3

test/MHD/Coarsening/CMakeLists.txt

@@ -0,0 +1 @@
+enable_idefix_property(Idefix_MHD)

test/MHD/Coarsening/definitions.hpp

@@ -0,0 +1,5 @@
+#define     COMPONENTS      3
+#define     DIMENSIONS      3
+
+
+#define     GEOMETRY        CARTESIAN

test/MHD/Coarsening/idefix-rkl.ini

@@ -0,0 +1,32 @@
+[Grid]
+X1-grid       1       -0.5  64  u  0.5
+X2-grid       1       -0.5  64  u  0.5
+X3-grid       1       -0.5  2   u  0.5
+coarsening    static  X1
+
+[TimeIntegrator]
+CFL         0.9
+tstop       0.1
+first_dt    1.e-4
+nstages     2
+
+[Hydro]
+solver         hlld
+resistivity    rkl   constant  0.05
+
+[Boundary]
+X1-beg    periodic
+X1-end    periodic
+X2-beg    periodic
+X2-end    periodic
+X3-beg    periodic
+X3-end    periodic
+
+[Output]
+vtk    0.1
+log    10
+# dmp         10.0
+
+[Setup]
+direction    1  0    # first is the field component, second is the direction of diffusion
+# advectionSpeed 0.2 # Advection speed

test/MHD/Coarsening/idefix-x2.ini

@@ -0,0 +1,32 @@
+[Grid]
+X1-grid       1       -0.5  64  u  0.5
+X2-grid       1       -0.5  64  u  0.5
+X3-grid       1       -0.5  2   u  0.5
+coarsening    static  X2
+
+[TimeIntegrator]
+CFL         0.9
+tstop       0.1
+first_dt    1.e-4
+nstages     2
+
+[Hydro]
+solver         hlld
+resistivity    explicit  constant  0.05
+
+[Boundary]
+X1-beg    periodic
+X1-end    periodic
+X2-beg    periodic
+X2-end    periodic
+X3-beg    periodic
+X3-end    periodic
+
+[Output]
+vtk    0.1
+log    10
+# dmp         10.0
+
+[Setup]
+direction    0  1    # first is the field component, second is the direction of diffusion
+# advectionSpeed 0.2 # Advection speed

test/MHD/Coarsening/idefix-x3.ini

@@ -0,0 +1,32 @@
+[Grid]
+X1-grid       1       -0.5  64  u  0.5
+X2-grid       1       -0.5  2   u  0.5
+X3-grid       1       -0.5  64  u  0.5
+coarsening    static  X3
+
+[TimeIntegrator]
+CFL         0.9
+tstop       0.1
+first_dt    1.e-4
+nstages     2
+
+[Hydro]
+solver         hlld
+resistivity    explicit  constant  0.05
+
+[Boundary]
+X1-beg    periodic
+X1-end    periodic
+X2-beg    periodic
+X2-end    periodic
+X3-beg    periodic
+X3-end    periodic
+
+[Output]
+vtk    0.1
+log    10
+# dmp         10.0
+
+[Setup]
+direction    0  2    # first is the field component, second is the direction of diffusion
+# advectionSpeed 0.2 # Advection speed

test/MHD/Coarsening/idefix.ini

@@ -0,0 +1,32 @@
+[Grid]
+X1-grid       1       -0.5  64  u  0.5
+X2-grid       1       -0.5  64  u  0.5
+X3-grid       1       -0.5  2   u  0.5
+coarsening    static  X1
+
+[TimeIntegrator]
+CFL         0.9
+tstop       0.1
+first_dt    1.e-4
+nstages     2
+
+[Hydro]
+solver         hlld
+resistivity    explicit  constant  0.05
+
+[Boundary]
+X1-beg    periodic
+X1-end    periodic
+X2-beg    periodic
+X2-end    periodic
+X3-beg    periodic
+X3-end    periodic
+
+[Output]
+vtk    0.1
+log    10
+dmp    0.1
+
+[Setup]
+direction    1  0    # first is the field component, second is the direction of diffusion
+# advectionSpeed 0.2 # Advection speed

test/MHD/Coarsening/python/testidefix.py

@@ -0,0 +1,96 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Mar  5 11:29:41 2020
+
+@author: glesur
+"""
+
+import os
+import sys
+sys.path.append(os.getenv("IDEFIX_DIR"))
+from pytools.vtk_io import readVTK
+import argparse
+import numpy as np
+import matplotlib.pyplot as plt
+import inifix
+
+parser = argparse.ArgumentParser()
+parser.add_argument("-noplot",
+                    default=False,
+                    help="disable plotting",
+                    action="store_true")
+
+
+args, unknown=parser.parse_known_args()
+
+# find which ini file was used
+with open('../idefix.0.log','r') as fp:
+  lines = fp.readlines()
+  for line in lines:
+    if line.find(".ini") != -1:
+      s=line.split()
+      inputFile=s[5][:-1]
+
+input = inifix.load('../'+inputFile)
+componentDir = input["Setup"]["direction"][0]
+spatialDir = input["Setup"]["direction"][1]
+
+# Convert componentDir
+componentName=["BX1","BX2","BX3"]
+
+
+V=readVTK('../data.0001.vtk', geometry='cartesian')
+
+# left_state and right_state set p, rho and u
+# geometry sets left boundary on 0., right boundary on 1 and initial
+# position of the shock xi on 0.5
+# t is the time evolution for which positions and states in tube should be calculated
+# gamma denotes specific heat
+# note that gamma and npts are default parameters (1.4 and 500) in solve function
+
+x=V.x
+if spatialDir==0:
+  qside=V.data[componentName[componentDir]][:,4,0]
+  qcenter=V.data[componentName[componentDir]][:,32,0]
+  x=V.x
+if spatialDir==1:
+  qside=V.data[componentName[componentDir]][4,:,0]
+  qcenter=V.data[componentName[componentDir]][32,:,0]
+  x=V.y
+if spatialDir==2:
+  qside=V.data[componentName[componentDir]][4,0,:]
+  qcenter=V.data[componentName[componentDir]][32,0,:]
+  x=V.z
+
+#initial condition in the code
+q0=0.*x+0.1*(np.abs(x)<0.1)
+
+# compute solution to diffusion equation using Fourier transforms
+qf0=np.fft.fft(q0)
+k2=(2.0*np.pi*q0.size*np.fft.fftfreq(q0.size))**2
+qf=qf0*np.exp(-k2*V.t*0.05)
+qtheory=np.real(np.fft.ifft(qf))
+
+if(not args.noplot):
+    plt.figure(1)
+    plt.plot(x,qside,'+',markersize=6,label='side')
+    plt.plot(x,qcenter,'+',markersize=6,label='center')
+
+    #plt.plot(x,q0)
+    plt.plot(x,qtheory,label='theory')
+    plt.legend()
+    plt.title(componentName[componentDir])
+
+
+
+    plt.ioff()
+    plt.show()
+
+errorCenter=np.mean((qcenter-qtheory)**2)
+errorSide=np.mean((qside-qtheory)**2)
+print("Error Side=%e, Error Center=%e"%(errorSide,errorCenter))
+assert(errorCenter<2e-5)
+assert(errorSide<2e-8)
+
+print("SUCCESS!")

test/MHD/Coarsening/setup.cpp

@@ -0,0 +1,94 @@
+#include "idefix.hpp"
+#include "setup.hpp"
+
+
+static int componentDir;     // Which field component varies
+static int spatialDir;  // along Which direction it varies
+static int coarseningDir; // coarsening direction (only one in this test!)
+static real advSpeed;
+
+void CoarsenFunction(DataBlock &data) {
+  IdefixArray2D<int> coarseningLevel = data.coarseningLevel[coarseningDir];
+  int dir = 0;
+  if(coarseningDir == componentDir) dir = spatialDir;
+  else if(coarseningDir == spatialDir) dir = componentDir;
+  IdefixArray1D<real> x = data.x[dir];
+
+  // Choose which index we should use for x, knowing that coarseningLevel lacks the direction coarseningDir
+  bool useSecond = (dir == KDIR);
+  if(dir == JDIR && coarseningDir == KDIR) useSecond = true;
+  idefix_for("set_coarsening", 0, coarseningLevel.extent(0), 0, coarseningLevel.extent(1),
+      KOKKOS_LAMBDA(int j,int i) {
+        int idx = i;
+        if(useSecond) idx = j;
+        coarseningLevel(j,i) = 1 + static_cast<int>(6*(0.5-fabs(x(idx))));
+        });
+}
+
+
+// Initialisation routine. Can be used to allocate
+// Arrays or variables which are used later on
+Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) {
+  componentDir = input.Get<int>("Setup", "direction", 0);
+  spatialDir = input.Get<int>("Setup", "direction", 1);
+  if(componentDir == spatialDir) {
+    IDEFIX_ERROR("The field component and the direction dependence cannot be identical");
+  }
+  advSpeed = input.GetOrSet<real>("Setup","advectionSpeed",0,0.0);
+  if(data.haveGridCoarsening) {
+    data.EnrollGridCoarseningLevels(&CoarsenFunction);
+    coarseningDir = -1;
+    int i = 0;
+    // detect coarsening direction
+    while(coarseningDir < 0) {
+      if(data.coarseningDirection[i]) coarseningDir = i;
+      i++;
+    }
+  }
+
+}
+
+// This routine initialize the flow
+// Note that data is on the device.
+// One can therefore define locally
+// a datahost and sync it, if needed
+void Setup::InitFlow(DataBlock &data) {
+    // Create a host copy
+    DataBlockHost d(data);
+
+    int idx=0;
+    for(int k = 0; k < d.np_tot[KDIR] ; k++) {
+      if(spatialDir == KDIR) idx = k;
+      for(int j = 0; j < d.np_tot[JDIR] ; j++) {
+        if(spatialDir == JDIR) idx = j;
+        for(int i = 0; i < d.np_tot[IDIR] ; i++) {
+          if(spatialDir == IDIR) idx = i;
+
+          real x=d.x[spatialDir](idx);
+
+          d.Vc(RHO,k,j,i) = 1.0;
+          d.Vc(VX1,k,j,i) = 0.0;
+          d.Vc(VX2,k,j,i) = 0.0;
+          d.Vc(VX3,k,j,i) = 0.0;
+
+          d.Vs(BX1s,k,j,i) = 0.0;
+          d.Vs(BX2s,k,j,i) = 0.0;
+          d.Vs(BX3s,k,j,i) = 0.0;
+
+          d.Vs(componentDir,k,j,i) = fabs(x) < 0.1 ? 0.1 : 0.0;
+
+          real Pmag = 0;
+          for(int n = 0 ; n < DIMENSIONS ; n++) {
+            Pmag += 0.5 * d.Vs(n,k,j,i) * d.Vs(n,k,j,i);
+          }
+
+          d.Vc(PRS,k,j,i) = 1.0 - Pmag;
+
+          d.Vc(VX1+spatialDir,k,j,i) = advSpeed;
+            }
+        }
+    }
+
+    // Send it all, if needed
+    d.SyncToDevice();
+}