Adiabatic expansion test ICs and documentation

docs/sphinx/ChollaExamples/3D-Adiabatic-Expansion/3D-Adiabatic-Expansion.md

@@ -0,0 +1,60 @@
+# 3D Adaibatic Expansion Test
+
+This test demonstrates simple adiabatic expansion in a cosmological frame.  The initial conditions are the critical density in baryons at the starting redshift, T=100K, and no bulk velocity. Gamma is set to 5/3. This test is performed with the cosmology build (`cholla/builds/make.type.cosmology`) and SIMPLE integrator.
+Full initial conditions can be found in `cholla/src/grid/initial_conditions.cpp`under `Adiabatic_Expansion()`. 
+
+## Parameter file:
+
+This parameter file can be found in [examples/3D/Adiabatic_Expansion.txt](https://github.com/cholla-hydro/cholla/blob/main/examples/3D/Adiabatic_Expansion.txt) on the `dev` branch.
+```
+#
+# Parameter File for the 3D Adiabatic Expansion test.
+#
+
+######################################
+# number of grid cells in the x dimension
+nx=256
+# number of grid cells in the y dimension
+ny=32
+# number of grid cells in the z dimension
+nz=32
+# output time
+tout=1000
+# how often to output
+outstep=1000
+# value of gamma
+gamma=1.66666667
+# name of initial conditions
+init=Adiabatic_Expansion
+#Cosmological Parameters
+Init_redshift=20.0
+#Init_redshift=0.998294693667
+H0=50.0
+Omega_M=1.0
+Omega_L=0.0
+Omega_b=1.0
+temperature_floor=1.0e-2
+scale_outputs_file=scale_output_files/outputs_adiabatic_expansion.txt
+# domain properties
+xmin=0.0
+ymin=0.0
+zmin=0.0
+xlen=64000.0
+ylen=8000.0
+zlen=8000.0
+# type of boundary conditions
+xl_bcnd=1
+xu_bcnd=1
+yl_bcnd=1
+yu_bcnd=1
+zl_bcnd=1
+zu_bcnd=1
+# path to output directory
+indir=ics/
+outdir=./
+```
+
+Upon completion, you should obtain two output files. The final density, velocity, and temperature in physical units are shown below. 
+:::{figure} adiabatic_expansion.png
+Adiabatic Expansion test solution from Cholla.
+::: 

docs/sphinx/ChollaExamples/index.md

@@ -9,5 +9,6 @@ Continue porting me from the [wiki](https://github.com/cholla-hydro/cholla/wiki/
 
 1D-123-Test/overview
 2D-Rayleigh-Taylor/2D-Rayleigh-Taylor
+3D-Adiabatic-Expansion/3D-Adiabatic-Expansion.md
 :::
 

examples/3D/Adiabatic_Expansion.txt

@@ -0,0 +1,45 @@
+#
+# Parameter File for the 3D Adiabatic Expansion test.
+#
+
+######################################
+# number of grid cells in the x dimension
+nx=256
+# number of grid cells in the y dimension
+ny=32
+# number of grid cells in the z dimension
+nz=32
+# output time
+tout=1000
+# how often to output
+outstep=1000
+# value of gamma
+gamma=1.66666667
+# name of initial conditions
+init=Adiabatic_Expansion
+#Cosmological Parameters 
+Init_redshift=20.0
+#Init_redshift=0.998294693667
+H0=50.0
+Omega_M=1.0
+Omega_L=0.0
+Omega_b=1.0
+temperature_floor=1.0e-2
+scale_outputs_file=scale_output_files/outputs_zeldovich_grav4.txt
+# domain properties
+xmin=0.0
+ymin=0.0
+zmin=0.0
+xlen=64000.0
+ylen=8000.0
+zlen=8000.0
+# type of boundary conditions
+xl_bcnd=1
+xu_bcnd=1
+yl_bcnd=1
+yu_bcnd=1
+zl_bcnd=1
+zu_bcnd=1
+# path to output directory
+indir=ics/
+outdir=./

src/grid/grid3D.h

@@ -694,6 +694,8 @@ class Grid3D
 
   void Zeldovich_Pancake(struct Parameters P);
 
+  void Adiabatic_Expansion(struct Parameters P);
+
   void Chemistry_Test(struct Parameters P);
 
 #ifdef MHD

src/grid/initial_conditions.cpp

@@ -82,6 +82,8 @@ void Grid3D::Set_Initial_Conditions(Parameters P)
     Uniform_Grid();
   } else if (strcmp(P.init, "Zeldovich_Pancake") == 0) {
     Zeldovich_Pancake(P);
+  } else if (strcmp(P.init, "Adiabatic_Expansion") == 0) {
+    Adiabatic_Expansion(P);
   } else if (strcmp(P.init, "Chemistry_Test") == 0) {
     Chemistry_Test(P);
 #ifdef MHD
@@ -1557,6 +1559,79 @@ void Grid3D::Zeldovich_Pancake(struct Parameters P)
         E    = ics_values[2 * nPoints + index];
         U    = ics_values[3 * nPoints + index];
         // //
+        // temp = T_init * pow( dens/rho_0, 2./3.);
+        // U = temp / (gamma-1) / MP * KB * 1e-10 * dens;
+        // E = 0.5*dens*vel*vel + U;
+
+        // chprintf( "%f \n", vel );
+        C.density[id]    = dens;
+        C.momentum_x[id] = dens * vel;
+        C.momentum_y[id] = 0;
+        C.momentum_z[id] = 0;
+        C.Energy[id]     = E;
+
+  #ifdef DE
+        C.GasEnergy[id] = U;
+  #endif
+      }
+    }
+  }
+
+#endif  // COSMOLOGY
+}
+
+void Grid3D::Adiabatic_Expansion(struct Parameters P)
+{
+#ifndef COSMOLOGY
+  chprintf("To run an Adiabatic Expansion test, COSMOLOGY has to be turned ON \n");
+  exit(-1);
+#else
+
+  int i, j, k, id;
+  Real x_pos, y_pos, z_pos;
+  Real H0, h, Omega_M, rho_0, G, z_zeldovich, z_init, x_center, T_init, k_x;
+
+  chprintf("Setting Adiabatic Expansion initial conditions...\n");
+  H0      = P.H0;
+  h       = H0 / 100;
+  Omega_M = P.Omega_M;
+
+  chprintf(" h = %f \n", h);
+  chprintf(" Omega_M = %f \n", Omega_M);
+
+  H0 /= 1000;  //[km/s / kpc]
+  G      = G_COSMO;
+  rho_0  = 3 * H0 * H0 / (8 * M_PI * G) * Omega_M / h / h;
+  z_init = P.Init_redshift;
+  chprintf(" rho_0 = %f \n", rho_0);
+  chprintf(" z_init = %f \n", z_init);
+
+  T_init = 100;
+  chprintf(" T initial = %f \n", T_init);
+
+  k_x = 2 * M_PI / H.xdglobal;
+
+  Real dens, vel, temp, U, E, gamma;
+  gamma = P.gamma;
+
+  int index;
+  // set the initial values of the conserved variables
+  for (k = H.n_ghost; k < H.nz - H.n_ghost; k++) {
+    for (j = H.n_ghost; j < H.ny - H.n_ghost; j++) {
+      for (i = H.n_ghost; i < H.nx - H.n_ghost; i++) {
+        id = i + j * H.nx + k * H.nx * H.ny;
+
+        // Analytical Initial Conditions
+        //  dens = rho_0 / ( 1 - ( 1 + z_zeldovich ) / ( 1 + z_init ) * cos(
+        //  k_x*( x_pos - x_center )) ); vel = - H0 * ( 1 + z_zeldovich ) /
+        //  sqrt( 1 + z_init ) * sin( k_x*( x_pos - x_center )) / k_x; temp =
+        //  T_init * pow( dens / rho_0, 2./3 ); U = temp / (gamma - 1) / MP * KB
+        //  * 1e-10 * dens; E = 0.5 * dens * vel * vel + U;
+
+        dens = rho_0;
+        vel  = 0;
+        U    = T_init / (gamma - 1) / MP * KB * 1e-10 * rho_0;
+        E    = U;
 
         // chprintf( "%f \n", vel );
         C.density[id]    = dens;

src/particles/particles_3D.cpp

@@ -214,6 +214,8 @@ void Particles3D::Initialize(Parameters *P, const SpatialDomainProps &spatial_pr
     Initialize_Sphere(P);
   } else if (strcmp(P->init, "Zeldovich_Pancake") == 0) {
     Initialize_Zeldovich_Pancake(P);
+  } else if (strcmp(P->init, "Adiabatic_Expansion") == 0) {
+    Initialize_Adiabatic_Expansion(P);
   } else if (strcmp(P->init, "Read_Grid") == 0) {
     Load_Particles_Data(P);
   } else if (strcmp(P->init, "Isolated_Stellar_Cluster") == 0) {
@@ -1139,6 +1141,17 @@ void Particles3D::Initialize_Zeldovich_Pancake(struct Parameters *P)
   chprintf(" Particles Zeldovich Pancake Initialized, n_local: %lu\n", n_local);
 }
 
+void Particles3D::Initialize_Adiabatic_Expansion(struct Parameters *P)
+{
+  // No particles for the Adiabatic Expansion problem. n_local=0
+
+  chprintf("Setting Adiabatic Expansion initial conditions...\n");
+
+  n_local = 0;
+
+  chprintf(" Particles Adiabatic Expansion initialized, n_local: %lu\n", n_local);
+}
+
 void Grid3D::Initialize_Uniform_Particles()
 {
   // Initialize positions assigning one particle at each cell in a uniform grid

src/particles/particles_3D.h

@@ -315,6 +315,8 @@ class Particles3D
 
   void Initialize_Zeldovich_Pancake(struct Parameters *P);
 
+  void Initialize_Adiabatic_Expansion(struct Parameters *P);
+
   void Load_Particles_Data(struct Parameters *P);
 
   void Free_Memory();