Added unit test for TRRS Riemann solver.

.gitignore

@@ -72,6 +72,7 @@ tests/testParser.sh
 tests/testReading.sh
 tests/testAdiabaticIndex
 tests/testRiemannExact
+tests/testRiemannTRRS
 
 theory/latex/swift.pdf
 theory/kernel/kernels.pdf

tests/Makefile.am

@@ -24,13 +24,13 @@ AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS
 TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetry \
         testPair.sh testPairPerturbed.sh test27cells.sh test27cellsPerturbed.sh  \
         testParser.sh testSPHStep test125cells.sh testKernelGrav testFFT \
-        testAdiabaticIndex testRiemannExact
+        testAdiabaticIndex testRiemannExact testRiemannTRRS
 
 # List of test programs to compile
 check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
 		 testSPHStep testPair test27cells test125cells testParser \
                  testKernel testKernelGrav testFFT testInteractions testMaths testSymmetry \
-                 testAdiabaticIndex testRiemannExact
+                 testAdiabaticIndex testRiemannExact testRiemannTRRS
 
 # Sources for the individual programs
 testGreetings_SOURCES = testGreetings.c
@@ -67,6 +67,8 @@ testAdiabaticIndex_SOURCES = testAdiabaticIndex.c
 
 testRiemannExact_SOURCES = testRiemannExact.c
 
+testRiemannTRRS_SOURCES = testRiemannTRRS.c
+
 # Files necessary for distribution
 EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
 	     test27cells.sh test27cellsPerturbed.sh tolerance.dat testParser.sh \

tests/testRiemannExact.c

@@ -212,7 +212,7 @@ void check_riemann_exact() {
 }
 
 /**
- * @brief Check adiabatic index constants and the various Riemann solvers
+ * @brief Check the exact Riemann solver
  */
 int main() {
 

tests/testRiemannTRRS.c

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include "error.h"
+#include "riemann/riemann_trrs.h"
+
+/**
+ * @brief Check that a and b are consistent (up to some error)
+ *
+ * @param a First value
+ * @param b Second value
+ * @param s String used to identify this check in messages
+ */
+void check_value(float a, float b, const char* s) {
+  if (fabsf(a - b) / fabsf(a + b) > 1.e-5f && fabsf(a - b) > 1.e-5f) {
+    error("Values are inconsistent: %g %g (%s)!", a, b, s);
+  } else {
+    message("Values are consistent: %g %g (%s).", a, b, s);
+  }
+}
+
+struct riemann_statevector {
+  /*! @brief Density */
+  float rho;
+
+  /*! @brief Fluid velocity */
+  float v;
+
+  /*! @brief Pressure */
+  float P;
+};
+
+/**
+ * @brief Check that the solution to the Riemann problem with given left and
+ * right state is consistent with the given expected solution
+ *
+ * @param WL Left state
+ * @param WR Right state
+ * @param Whalf Expected solution
+ * @param s String used to identify this check in messages
+ */
+void check_riemann_solution(struct riemann_statevector* WL,
+                            struct riemann_statevector* WR,
+                            struct riemann_statevector* Whalf, const char* s) {
+  float WLarr[5], WRarr[5], Whalfarr[5], n_unit[3];
+
+  n_unit[0] = 1.0f;
+  n_unit[1] = 0.0f;
+  n_unit[2] = 0.0f;
+
+  WLarr[0] = WL->rho;
+  WLarr[1] = WL->v;
+  WLarr[2] = 0.0f;
+  WLarr[3] = 0.0f;
+  WLarr[4] = WL->P;
+
+  WRarr[0] = WR->rho;
+  WRarr[1] = WR->v;
+  WRarr[2] = 0.0f;
+  WRarr[3] = 0.0f;
+  WRarr[4] = WR->P;
+
+  riemann_solver_solve(WLarr, WRarr, Whalfarr, n_unit);
+
+  message("Checking %s...", s);
+  check_value(Whalfarr[0], Whalf->rho, "rho");
+  check_value(Whalfarr[1], Whalf->v, "v");
+  check_value(Whalfarr[4], Whalf->P, "P");
+}
+
+/**
+ * @brief Check the TRRS Riemann solver on the Toro test problems
+ */
+void check_riemann_trrs() {
+  struct riemann_statevector WL, WR, Whalf;
+
+  /* Test 1 */
+  WL.rho = 1.0f;
+  WL.v = 0.0f;
+  WL.P = 1.0f;
+  WR.rho = 0.125f;
+  WR.v = 0.0f;
+  WR.P = 0.1f;
+#if defined(HYDRO_GAMMA_5_3)
+  Whalf.rho = 0.481167f;
+  Whalf.v = 0.838085f;
+  Whalf.P = 0.295456f;
+#elif defined(HYDRO_GAMMA_4_3)
+  Whalf.rho = 0.41586f;
+  Whalf.v = 0.942546f;
+  Whalf.P = 0.310406f;
+#elif defined(HYDRO_GAMMA_2_1)
+  Whalf.rho = 0.53478f;
+  Whalf.v = 0.760037f;
+  Whalf.P = 0.285989f;
+#else
+#error "Unsupported adiabatic index!"
+#endif
+  check_riemann_solution(&WL, &WR, &Whalf, "Test 1");
+
+  /* Test 2 */
+  WL.rho = 1.0f;
+  WL.v = -2.0f;
+  WL.P = 0.4f;
+  WR.rho = 1.0f;
+  WR.v = 2.0f;
+  WR.P = 0.4f;
+#if defined(HYDRO_GAMMA_5_3)
+  Whalf.rho = 0.00617903f;
+  Whalf.v = 0.0f;
+  Whalf.P = 8.32249e-5f;
+#elif defined(HYDRO_GAMMA_4_3)
+  Whalf.rho = 0.0257933f;
+  Whalf.v = 0.0f;
+  Whalf.P = 0.00304838f;
+#elif defined(HYDRO_GAMMA_2_1)
+  Whalf.rho = 0.013932f;
+  Whalf.v = 0.0f;
+  Whalf.P = 7.76405e-5f;
+#else
+#error "Unsupported adiabatic index!"
+#endif
+  check_riemann_solution(&WL, &WR, &Whalf, "Test 2");
+
+  /* Test 3 */
+  WL.rho = 1.0f;
+  WL.v = 0.0f;
+  WL.P = 1000.0f;
+  WR.rho = 1.0f;
+  WR.v = 0.0f;
+  WR.P = 0.01f;
+#if defined(HYDRO_GAMMA_5_3)
+  Whalf.rho = 0.919498f;
+  Whalf.v = 3.37884f;
+  Whalf.P = 869.464f;
+#elif defined(HYDRO_GAMMA_4_3)
+  Whalf.rho = 0.941258f;
+  Whalf.v = 2.19945f;
+  Whalf.P = 922.454f;
+#elif defined(HYDRO_GAMMA_2_1)
+  Whalf.rho = 0.902032f;
+  Whalf.v = 4.49417f;
+  Whalf.P = 813.662f;
+#else
+#error "Unsupported adiabatic index!"
+#endif
+  check_riemann_solution(&WL, &WR, &Whalf, "Test 3");
+
+  /* Test 4 */
+  WL.rho = 1.0f;
+  WL.v = 0.0f;
+  WL.P = 0.01f;
+  WR.rho = 1.0f;
+  WR.v = 0.0f;
+  WR.P = 100.0f;
+#if defined(HYDRO_GAMMA_5_3)
+  Whalf.rho = 0.857525f;
+  Whalf.v = -1.93434f;
+  Whalf.P = 77.4007f;
+#elif defined(HYDRO_GAMMA_4_3)
+  Whalf.rho = 0.880649f;
+  Whalf.v = -1.45215f;
+  Whalf.P = 84.4119f;
+#elif defined(HYDRO_GAMMA_2_1)
+  Whalf.rho = 0.843058f;
+  Whalf.v = -2.31417f;
+  Whalf.P = 71.0747f;
+#else
+#error "Unsupported adiabatic index!"
+#endif
+  check_riemann_solution(&WL, &WR, &Whalf, "Test 4");
+
+  /* Test 5 */
+  WL.rho = 5.99924f;
+  WL.v = 19.5975f;
+  WL.P = 460.894f;
+  WR.rho = 5.99242f;
+  WR.v = -6.19633f;
+  WR.P = 46.0950f;
+#if defined(HYDRO_GAMMA_5_3)
+  Whalf.rho = 5.99924f;
+  Whalf.v = 19.5975f;
+  Whalf.P = 460.894f;
+#elif defined(HYDRO_GAMMA_4_3)
+  Whalf.rho = 5.99924f;
+  Whalf.v = 19.5975f;
+  Whalf.P = 460.894f;
+#elif defined(HYDRO_GAMMA_2_1)
+  Whalf.rho = 5.99924f;
+  Whalf.v = 19.5975f;
+  Whalf.P = 460.894f;
+#else
+#error "Unsupported adiabatic index!"
+#endif
+  check_riemann_solution(&WL, &WR, &Whalf, "Test 5");
+}
+
+/**
+ * @brief Check the TRRS Riemann solver
+ */
+int main() {
+
+  /* check the TRRS Riemann solver */
+  check_riemann_trrs();
+
+  return 0;
+}