add 1D reference case, and new setup

Author: glesur

test/HD/SedovBlastWave/idefix.ini

@@ -1,18 +1,21 @@
 [Grid]
-X1-grid    1  -0.5  256  u  0.5
-X2-grid    1  -0.5  256  u  0.5
-X3-grid    1  -0.5  256  u  0.5
+X1-grid    1  -0.5  128  u  0.5
+X2-grid    1  -0.5  128  u  0.5
+X3-grid    1  -0.5  128  u  0.5
 
 [TimeIntegrator]
 CFL         0.9
-tstop       5.0
+tstop       0.1
 first_dt    1.e-6
 nstages     2
 
 [Hydro]
 solver    hllc
 gamma     1.666666666666666666
 
+[Setup]
+Rstart    0.03
+
 [Boundary]
 X1-beg    periodic
 X1-end    periodic
@@ -22,4 +25,4 @@ X3-beg    periodic
 X3-end    periodic
 
 [Output]
-vtk    0.01
+vtk    0.1

test/HD/SedovBlastWave/python/1Dsolution/definitions.hpp

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

test/HD/SedovBlastWave/python/1Dsolution/idefix.ini

@@ -0,0 +1,26 @@
+[Grid]
+X1-grid    1  0.001  4096  u  0.5
+X2-grid    1  -0.5   1     u  0.5
+X3-grid    1  -0.5   1     u  0.5
+
+[TimeIntegrator]
+CFL         0.9
+tstop       0.1
+first_dt    1.e-9
+nstages     2
+
+[Hydro]
+solver    hllc
+gamma     1.666666666666666666
+
+[Boundary]
+X1-beg    outflow
+X1-end    outflow
+X2-beg    periodic
+X2-end    periodic
+X3-beg    periodic
+X3-end    periodic
+
+[Output]
+vtk    0.1
+log    1000

test/HD/SedovBlastWave/python/1Dsolution/setup.cpp

@@ -0,0 +1,52 @@
+#include "idefix.hpp"
+#include "setup.hpp"
+
+
+// Default constructor
+
+
+// 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) {
+
+}
+
+// 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++) {
+              real r = d.x[IDIR](i);
+
+                // Sedov Blast Wave Following Stone+2018, 3.4.4
+                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.Vc(PRS,k,j,i) = 0.01;
+
+                if(r<0.01) {
+                  d.Vc(PRS,k,j,i) = 1.6e5;
+                }
+
+
+            }
+        }
+    }
+
+    // Send it all, if needed
+    d.SyncToDevice();
+}
+
+// Analyse data to produce an output
+void MakeAnalysis(DataBlock & data) {
+
+}

test/HD/SedovBlastWave/python/testidefix.py

@@ -0,0 +1,61 @@
+#!/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
+from pytools import sod
+import argparse
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.interpolate import interp1d
+
+parser = argparse.ArgumentParser()
+parser.add_argument("-noplot",
+                    default=False,
+                    help="disable plotting",
+                    action="store_true")
+
+
+args, unknown=parser.parse_known_args()
+
+V=readVTK('../data.0001.vtk')
+R=readVTK('1Dsolution/data.0001.vtk')
+
+
+
+# Finally, let's plot solutions
+p = R.data['PRS'][:,0,0]
+x= R.r
+
+[x3D, y3D, z3D]= np.meshgrid(V.x,V.y,V.z,indexing='ij')
+r3D=np.sqrt(x3D**2+y3D**2+z3D**2)
+p3D=V.data['PRS']
+
+
+solinterp=interp1d(x,p)
+
+
+if(not args.noplot):
+    plt.figure(1)
+    plt.plot(x,p)
+    plt.plot(r3D.flatten(),p3D.flatten(),'+',markersize=2)
+    plt.title('Pressure')
+
+    plt.ioff()
+    plt.show()
+
+error=np.mean(np.fabs(V.data['PRS'][:,0,0]-solinterp(V.x)))
+print("Error=%e"%error)
+if error<2e-3:
+    print("SUCCESS!")
+    sys.exit(0)
+else:
+    print("FAILURE!")
+    sys.exit(1)

test/HD/SedovBlastWave/setup.cpp

@@ -4,11 +4,11 @@
 
 // Default constructor
 
-
+real Rstart;
 // 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) {
-
+  Rstart = input.Get<real>("Setup","Rstart",0);
 }
 
 // This routine initialize the flow
@@ -19,26 +19,44 @@ void Setup::InitFlow(DataBlock &data) {
     // Create a host copy
     DataBlockHost d(data);
 
+    real V = 0;
+    for(int k = 0; k < d.np_int[KDIR] ; k++) {
+        for(int j = 0; j < d.np_int[JDIR] ; j++) {
+            for(int i = 0; i < d.np_int[IDIR] ; i++) {
+              real x = d.x[IDIR](i);
+              real y = d.x[JDIR](j);
+              real z = d.x[KDIR](k);
+              real r=sqrt(x*x+y*y+z*z);
+
+              if(r<Rstart) {
+                V += d.dV(k,j,i);
+              }
+    }}}
+
+#ifdef WITH_MPI
+    MPI_Allreduce(MPI_IN_PLACE, &V, 1, realMPI, MPI_SUM, MPI_COMM_WORLD);
+#endif
 
+    real gamma = data.hydro.GetGamma();
     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++) {
               real x = d.x[IDIR](i);
               real y = d.x[JDIR](j);
               real z = d.x[KDIR](k);
 
-                // Sedov Blast Wave Following Stone+2018, 3.4.4
-                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;
+              // Sedov Blast Wave Following Stone+2018, 3.4.4
+              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.Vc(PRS,k,j,i) = 0.01;
+              d.Vc(PRS,k,j,i) = 0.01;
 
-                real r=sqrt(x*x+y*y+z*z);
-                if(r<0.01) {
-                  d.Vc(PRS,k,j,i) = 1.6e5;
-                }
+              real r=sqrt(x*x+y*y+z*z);
+              if(r<Rstart) {
+                d.Vc(PRS,k,j,i) = (gamma-1)/V;
+              }
 
 
             }