regression test for thermal diffusion

Author: glesur

test/HD/thermalDiffusion/definitions.hpp

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

test/HD/thermalDiffusion/idefix-rkl.ini

@@ -0,0 +1,35 @@
+[Grid]
+
+X1-grid    1    -0.5    500    u    0.5
+X2-grid    1    0.0    1    u    1.0
+X3-grid    1    0.0    1      u    1.0
+
+
+
+[TimeIntegrator]
+
+CFL              0.8
+CFL_max_var      1.1
+tstop            0.2
+nstages          2
+
+[Hydro]
+
+solver         hllc
+gamma          1.4
+TDiffusion  rkl  constant  0.1
+
+[Setup]
+amplitude     1e-6
+
+[Boundary]
+
+X1-beg        periodic
+X1-end        periodic
+X2-beg        periodic
+X2-end        periodic
+X3-beg        periodic
+X3-end        periodic
+
+[Output]
+analysis    0.01

test/HD/thermalDiffusion/idefix.ini

@@ -0,0 +1,35 @@
+[Grid]
+
+X1-grid    1    -0.5    500    u    0.5
+X2-grid    1    0.0    1    u    1.0
+X3-grid    1    0.0    1      u    1.0
+
+
+
+[TimeIntegrator]
+
+CFL              0.8
+CFL_max_var      1.1
+tstop            0.2
+nstages          2
+
+[Hydro]
+
+solver         hllc
+gamma          1.4
+TDiffusion  explicit  constant  0.1
+
+[Setup]
+amplitude     1e-6
+
+[Boundary]
+
+X1-beg        periodic
+X1-end        periodic
+X2-beg        periodic
+X2-end        periodic
+X3-beg        periodic
+X3-end        periodic
+
+[Output]
+analysis    0.01

test/HD/thermalDiffusion/python/testidefix.py

@@ -0,0 +1,52 @@
+import sys
+import numpy as np
+import argparse
+import matplotlib.pyplot as plt
+from scipy.signal import argrelextrema
+
+parser = argparse.ArgumentParser()
+parser.add_argument("-noplot",
+                    default=False,
+                    help="disable plotting",
+                    action="store_true")
+
+
+args=parser.parse_args()
+
+
+# values from the setup to compute the decay rate
+kappa=0.1
+gamma=1.4
+rho0=1.0
+
+#diffusion coefficient
+D=kappa*(gamma-1)/rho0
+
+k=2.0*np.pi
+
+# decay rate of the temperature perturbation
+gamma=-k**2*D
+
+raw=np.loadtxt('../analysis.dat',skiprows=1)
+t=raw[:,0]
+Tidfx=raw[:,1]
+
+Tth=Tidfx[0]*np.exp(gamma*t)
+error=np.mean(Tidfx/Tth-1.0)
+
+if(not args.noplot):
+
+    plt.close('all')
+    plt.figure()
+    plt.plot(t,Tidfx/Tth-1.0)
+
+    print("Error=%e"%error)
+    plt.ioff()
+    plt.show()
+
+if(error<2.5e-6):
+    print("SUCCESS")
+    sys.exit(0)
+else:
+    print("Failed")
+    sys.exit(1)

test/HD/thermalDiffusion/setup.cpp

@@ -0,0 +1,80 @@
+#include "idefix.hpp"
+#include "setup.hpp"
+
+
+#define FILENAME  "analysis.dat"
+real amplitude;
+// Analyse data to get the amplitude of the k=1 temperature perturbation
+void Analysis(DataBlock & data) {
+  DataBlockHost d(data);
+  d.SyncFromDevice();
+  real Tmode = 0;
+  for(int i = d.beg[IDIR]; i < d.end[IDIR] ; i++) {
+    real T = d.Vc(PRS,d.beg[KDIR],d.beg[JDIR],i) / d.Vc(RHO,d.beg[KDIR],d.beg[JDIR],i);
+    Tmode += T * sin(2.0*M_PI*d.x[IDIR](i));
+  }
+  Tmode = 2*Tmode/d.np_int[JDIR];
+
+  std::ofstream f;
+  f.open(FILENAME,std::ios::app);
+  f.precision(10);
+  f << std::scientific << data.t << "\t" << Tmode << std::endl;
+  f.close();
+
+}
+
+void InternalBoundary(DataBlock& data, const real t) {
+  IdefixArray4D<real> Vc = data.hydro.Vc;
+  idefix_for("InternalBoundary",0,data.np_tot[KDIR],0,data.np_tot[JDIR],0,data.np_tot[IDIR],
+              KOKKOS_LAMBDA (int k, int j, int i) {
+                // Cancel any motion that could be happening
+                Vc(VX1,k,j,i) = 0.0;
+                Vc(VX2,k,j,i) = 0.0;
+                Vc(VX3,k,j,i) = 0.0;
+              });
+}
+
+
+// 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) {
+  output.EnrollAnalysis(&Analysis);
+  data.hydro.EnrollInternalBoundary(&InternalBoundary);
+
+  amplitude = input.GetReal("Setup","amplitude",0);
+  // Initialise the output file
+  std::ofstream f;
+  f.open(FILENAME,std::ios::trunc);
+  f << "t\t\t T" << std::endl;
+  f.close();
+}
+
+// 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);
+
+
+    for(int k = 0; k < d.np_tot[KDIR] ; k++) {
+        for(int j = 0; j < d.np_tot[JDIR] ; j++) {
+            for(int i = 0; i < d.np_tot[IDIR] ; i++) {
+
+                d.Vc(RHO,k,j,i) = 1 - amplitude*sin(2.0*M_PI*d.x[IDIR](i));
+                d.Vc(VX1,k,j,i) = ZERO_F;
+                d.Vc(PRS,k,j,i) = 1.0;
+
+            }
+        }
+    }
+
+    // Send it all, if needed
+    d.SyncToDevice();
+}
+
+// Analyse data to produce an output
+void MakeAnalysis(DataBlock & data) {
+
+}