New test to verify Euler source term

Checks source term against a finite-difference of the convective fluxes. Since the source term is analytically derived for the manufactured solution, it also checks for the physical flux in a way. Also, new test for grid convegence of Euler. Not supposed to pass Test project /home/ddong/Codes/PHiLiP/build_release Start 1: 1D_numerical_flux_conservation 1/27 Test #1: 1D_numerical_flux_conservation .............................. Passed 0.17 sec Start 2: 2D_numerical_flux_conservation 2/27 Test #2: 2D_numerical_flux_conservation .............................. Passed 0.17 sec Start 3: 3D_numerical_flux_conservation 3/27 Test #3: 3D_numerical_flux_conservation .............................. Passed 0.16 sec Start 4: 1D_jacobian_matrix_regression 4/27 Test #4: 1D_jacobian_matrix_regression ............................... Passed 0.45 sec Start 5: 2D_jacobian_matrix_regression 5/27 Test #5: 2D_jacobian_matrix_regression ............................... Passed 0.50 sec Start 6: 3D_jacobian_matrix_regression 6/27 Test #6: 3D_jacobian_matrix_regression ............................... Passed 1.57 sec Start 7: 1D_euler_convert_primitive_conservative 7/27 Test #7: 1D_euler_convert_primitive_conservative ..................... Passed 0.16 sec Start 8: 2D_euler_convert_primitive_conservative 8/27 Test #8: 2D_euler_convert_primitive_conservative ..................... Passed 0.20 sec Start 9: 3D_euler_convert_primitive_conservative 9/27 Test #9: 3D_euler_convert_primitive_conservative ..................... Passed 0.44 sec Start 10: 1D_euler_manufactured_solution_source 10/27 Test #10: 1D_euler_manufactured_solution_source ....................... Passed 0.16 sec Start 11: 2D_euler_manufactured_solution_source 11/27 Test #11: 2D_euler_manufactured_solution_source ....................... Passed 0.21 sec Start 12: 3D_euler_manufactured_solution_source 12/27 Test #12: 3D_euler_manufactured_solution_source ....................... Passed 0.54 sec Start 13: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION 13/27 Test #13: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ................. Passed 0.49 sec Start 14: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION 14/27 Test #14: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ................. Passed 2.38 sec Start 15: 3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION 15/27 Test #15: 3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ................. Passed 25.05 sec Start 16: 1D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION 16/27 Test #16: 1D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ................. Passed 0.58 sec Start 17: 2D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION 17/27 Test #17: 2D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ................. Passed 1.49 sec Start 18: 3D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION 18/27 Test #18: 3D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ................. Passed 48.04 sec Start 19: 1D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION 19/27 Test #19: 1D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ...... Passed 0.51 sec Start 20: 2D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION 20/27 Test #20: 2D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ...... Passed 1.50 sec Start 21: 3D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION 21/27 Test #21: 3D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ...... Passed 43.90 sec Start 22: 1D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION 22/27 Test #22: 1D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION ... Passed 0.49 sec Start 23: 2D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION 23/27 Test #23: 2D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION ... Passed 6.13 sec Start 24: 3D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION 24/27 Test #24: 3D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION ... Passed 20.49 sec Start 25: 1D_EULER_IMPLICIT_MANUFACTURED_SOLUTION 25/27 Test #25: 1D_EULER_IMPLICIT_MANUFACTURED_SOLUTION .....................***Failed 0.42 sec Start 26: 2D_EULER_IMPLICIT_MANUFACTURED_SOLUTION 26/27 Test #26: 2D_EULER_IMPLICIT_MANUFACTURED_SOLUTION .....................***Failed 0.41 sec Start 27: 3D_EULER_IMPLICIT_MANUFACTURED_SOLUTION 27/27 Test #27: 3D_EULER_IMPLICIT_MANUFACTURED_SOLUTION .....................***Failed 0.41 sec 89% tests passed, 3 tests failed out of 27 Total Test time (real) = 157.22 sec The following tests FAILED: 25 - 1D_EULER_IMPLICIT_MANUFACTURED_SOLUTION (Failed) 26 - 2D_EULER_IMPLICIT_MANUFACTURED_SOLUTION (Failed) 27 - 3D_EULER_IMPLICIT_MANUFACTURED_SOLUTION (Failed) Errors while running CTest
dougshidong · Jun 4, 2019 · 011a7c5 · 011a7c5
1 parent 5ea8857
commit 011a7c5
Show file tree

Hide file tree

Showing 9 changed files with 372 additions and 31 deletions.
diff --git a/src/physics/euler.cpp b/src/physics/euler.cpp
@@ -21,16 +21,17 @@ ::manufactured_solution (const dealii::Point<dim,double> &pos) const
     const double pi = atan(1)*4.0;
     const double ee = exp(1);
 
-    dealii::Table<1,real> base_value(nstate);
-    base_value[0] = pi/4;
+    std::array<real,nstate> base_value;
+    base_value[0] = pi/4.0;
     for (int i=0;i<dim;i++) {
         base_value[1+i] = ee/(i+1);
     }
     base_value[nstate-1] = ee/pi;
 
 
-    dealii::Table<2,real> amplitudes(nstate, dim);
-    dealii::Table<2,real> frequencies(nstate, dim);
+
+    std::array<dealii::Tensor<1,dim,real>,nstate> amplitudes;
+    std::array<dealii::Tensor<1,dim,real>,nstate> frequencies;
     for (int s=0; s<nstate; s++) {
         for (int d=0; d<dim; d++) {
             amplitudes[s][d] = 0.5*base_value[s]*(dim-d)/dim*(nstate-s)/nstate;
@@ -59,23 +60,85 @@ ::source_term (
     const double pi = atan(1)*4.0;
     const double ee = exp(1);
 
-    dealii::Table<1,real> base_value(nstate);
-    base_value[0] = pi/4;
+    std::array<real,nstate> base_value;
+    base_value[0] = pi/4.0;
     for (int i=0;i<dim;i++) {
         base_value[1+i] = ee/(i+1);
     }
     base_value[nstate-1] = ee/pi;
 
-    dealii::Table<2,real> amplitudes(nstate, dim);
-    dealii::Table<2,real> frequencies(nstate, dim);
+
+    std::array<dealii::Tensor<1,dim,real>,nstate> ampl;
+    std::array<dealii::Tensor<1,dim,real>,nstate> freq;
     for (int s=0; s<nstate; s++) {
-        for (int d=0; d<nstate; d++) {
-            amplitudes[s][d] = 0.5*base_value[s]*(dim-d)/dim*(nstate-s)/nstate;
-            frequencies[s][d] = 1.0+sin(0.1+s*0.5+d*0.2);
+        for (int d=0; d<dim; d++) {
+            ampl[s][d] = 0.5*base_value[s]*(dim-d)/dim*(nstate-s)/nstate;
+            freq[s][d] = 1.0+sin(0.1+s*0.5+d*0.2);
         }
     }
 
+
+    // Obtain primitive solution
+    const std::array<real,nstate> conservative_manu_soln = manufactured_solution(pos);
+    const std::array<real,nstate> primitive_soln = Euler<dim,nstate,real>::convert_conservative_to_primitive(conservative_manu_soln);
+    const real density = primitive_soln[0];
+    const std::array<real,dim> vel = Euler::extract_velocities_from_primitive(primitive_soln);
+    const real pressure = primitive_soln[nstate-1];
+    const real energy = conservative_manu_soln[nstate-1];
+
+    // Derivatives of total energy w.r.t primitive variables
+    const real v2 = Euler::compute_velocity_squared(vel);
+    const real dedr = 0.5*v2;
+    std::array<real,dim> dedv;
+    for (int d=0;d<dim;d++) {
+        dedv[d] = density*vel[d];
+    }
+    const real dedp = 1.0/(gam-1.0);
+
     std::array<real,nstate> source;
+    source.fill(0.0);
+    dealii::Table<3,real> dflux_dx(nstate, dim, dim); // state, flux_dim, deri_dim
+    for (int deri_dim=0; deri_dim<dim; ++deri_dim){
+
+        // Derivative of primitive solution w.r.t x, y, or z
+        real drdx = ampl[0][deri_dim] * freq[0][deri_dim] * pi/2.0 * cos( freq[0][deri_dim] * pos[deri_dim] * pi / 2.0 );
+        std::array<real,dim> dvdx;
+        for (int vel_d=0;vel_d<dim;vel_d++) {
+            dvdx[vel_d] = ampl[1+vel_d][deri_dim] * freq[1+vel_d][deri_dim] * pi/2.0 * cos( freq[1+vel_d][deri_dim] * pos[deri_dim] * pi / 2.0 );
+        }
+        real dpdx = ampl[nstate-1][deri_dim] * freq[nstate-1][deri_dim] * pi/2.0 * cos( freq[nstate-1][deri_dim] * pos[deri_dim] * pi / 2.0 );
+
+
+        //for (int flux_dim=0; flux_dim<dim; ++flux_dim){
+        // Only need the divergence of the flux, not all the derivatives
+        const int flux_dim = deri_dim;
+        {
+            const real vflux = primitive_soln[1+flux_dim];
+            const real dvflux_dx = dvdx[flux_dim];
+
+            const real dmassflux_dx = drdx * vflux + dvflux_dx * density;
+
+            // Mass flux
+            dflux_dx[0][flux_dim][deri_dim] = dmassflux_dx;
+            // Momentum flux
+            for (int vel_d=0; vel_d<dim; ++vel_d) {
+                dflux_dx[1+vel_d][flux_dim][deri_dim] = dvdx[vel_d] * density * vflux + dmassflux_dx*vel[vel_d];
+            }
+            dflux_dx[1+flux_dim][flux_dim][deri_dim] += dpdx;
+            // Energy flux
+            dflux_dx[nstate-1][flux_dim][deri_dim] = dvflux_dx * (energy + pressure) + dpdx * vflux;
+            dflux_dx[nstate-1][flux_dim][deri_dim] += dedr * drdx * vflux;
+            for (int vel_d=0; vel_d<dim; ++vel_d) {
+                dflux_dx[nstate-1][flux_dim][deri_dim] += dedv[vel_d] * dvdx[vel_d] * vflux;
+            }
+            dflux_dx[nstate-1][flux_dim][deri_dim] += dedp * dpdx * vflux;
+        }
+        source[0] += dflux_dx[0][flux_dim][deri_dim];
+        for (int vel_d=0; vel_d<dim; ++vel_d) {
+            source[1+vel_d] += dflux_dx[1+vel_d][flux_dim][deri_dim];
+        }
+        source[nstate-1] += dflux_dx[nstate-1][flux_dim][deri_dim];
+    }
 
     return source;
 }
@@ -90,11 +153,11 @@ ::convert_conservative_to_primitive ( const std::array<real,nstate> &conservativ
     std::array<real, dim> vel = compute_velocities (conservative_soln);
     real pressure = compute_pressure (conservative_soln);
 
-    primitive_soln[0] = conservative_soln[0];
+    primitive_soln[0] = density;
     for (int d=0; d<dim; ++d) {
         primitive_soln[1+d] = vel[d];
     }
-    primitive_soln[1+dim] = pressure;
+    primitive_soln[nstate-1] = pressure;
     return primitive_soln;
 }
 
@@ -105,14 +168,13 @@ ::convert_primitive_to_conservative ( const std::array<real,nstate> &primitive_s
 
     const real density = primitive_soln[0];
     const std::array<real,dim> velocities = extract_velocities_from_primitive(primitive_soln);
-    const real pressure = primitive_soln[1+dim];
 
     std::array<real, nstate> conservative_soln;
     conservative_soln[0] = density;
     for (int d=0; d<dim; ++d) {
         conservative_soln[1+d] = density*velocities[d];
     }
-    conservative_soln[1+dim] = compute_energy(primitive_soln);
+    conservative_soln[nstate-1] = compute_energy(primitive_soln);
 
     return conservative_soln;
 }
@@ -152,7 +214,7 @@ inline real Euler<dim,nstate,real>
 ::compute_energy ( const std::array<real,nstate> &primitive_soln ) const
 {
     const real density = primitive_soln[0];
-    const real pressure = primitive_soln[1+dim];
+    const real pressure = primitive_soln[nstate-1];
     const std::array<real,dim> velocities = extract_velocities_from_primitive(primitive_soln);
     const real vel2 = compute_velocity_squared(velocities);
 
@@ -165,7 +227,7 @@ inline real Euler<dim,nstate,real>
 ::compute_pressure ( const std::array<real,nstate> &conservative_soln ) const
 {
     const real density = conservative_soln[0];
-    const real energy  = conservative_soln[1+dim];
+    const real energy  = conservative_soln[nstate-1];
     const std::array<real,dim> vel = compute_velocities(conservative_soln);
     const real vel2 = compute_velocity_squared(vel);
     const real pressure = (gam-1.0)*(energy - 0.5*density*vel2);
@@ -192,16 +254,16 @@ ::convective_flux (const std::array<real,nstate> &conservative_soln) const
     const std::array<real,dim> vel = compute_velocities(conservative_soln);
     const real tot_energy = conservative_soln[nstate-1];
 
-    for (int fdim=0; fdim<dim; ++fdim) {
+    for (int flux_dim=0; flux_dim<dim; ++flux_dim) {
         // Density equation
-        conv_flux[0][fdim] = conservative_soln[1+fdim];
+        conv_flux[0][flux_dim] = conservative_soln[1+flux_dim];
         // Momentum equation
-        for (int sdim=0; sdim<dim; ++sdim){
-            conv_flux[1+sdim][fdim] = density*vel[fdim]*vel[sdim];
+        for (int velocity_dim=0; velocity_dim<dim; ++velocity_dim){
+            conv_flux[1+velocity_dim][flux_dim] = density*vel[flux_dim]*vel[velocity_dim];
         }
-        conv_flux[1+fdim][1+fdim] += pressure; // Add diagonal of pressure
+        conv_flux[1+flux_dim][flux_dim] += pressure; // Add diagonal of pressure
         // Energy equation
-        conv_flux[2+dim][fdim] = (conservative_soln[2+dim]+pressure)*vel[fdim];
+        conv_flux[nstate-1][flux_dim] = (tot_energy+pressure)*vel[flux_dim];
     }
     return conv_flux;
 }
@@ -214,6 +276,8 @@ ::convective_eigenvalues (
 {
     const std::array<real,dim> vel = compute_velocities(conservative_soln);
     std::array<real,nstate> eig;
+    (void) vel;
+    (void) normal;
     for (int i=0; i<nstate; i++) {
         //eig[i] = advection_speed*normal;
     }

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
@@ -1,11 +1,12 @@
+# Unit tests
 add_subdirectory(numerical_flux)
 add_subdirectory(regression)
+add_subdirectory(euler_unit_test)
 
+# Integration tests
 add_subdirectory(advection_implicit)
 add_subdirectory(diffusion_implicit)
 add_subdirectory(convection_diffusion_implicit)
-
 add_subdirectory(advection_vector_implicit)
 
-
-add_subdirectory(euler_unit_test)
+add_subdirectory(euler_implicit)
diff --git a/tests/euler_implicit/1d_euler_implicit.prm b/tests/euler_implicit/1d_euler_implicit.prm
@@ -0,0 +1,43 @@
+# Listing of Parameters
+# ---------------------
+# Number of dimensions
+set dimension = 1
+
+# The PDE we want to solve. Choices are
+# <diffusion|diffusion|convection_diffusion>.
+set pde_type  = euler
+
+subsection ODE solver
+
+  set ode_output                          = verbose
+
+  # Maximum nonlinear solver iterations
+  set nonlinear_max_iterations            = 500000
+
+  # Nonlinear solver residual tolerance
+  set nonlinear_steady_residual_tolerance = 1e-12
+
+  # Print every print_iteration_modulo iterations of the nonlinear solver
+  set print_iteration_modulo              = 1
+
+  # Explicit or implicit solverChoices are <explicit|implicit>.
+  set ode_solver_type                         = implicit
+end
+
+subsection manufactured solution convergence study
+  # Last degree used for convergence study
+  set degree_end        = 3
+
+  # Starting degree for convergence study
+  set degree_start      = 1
+
+  # Multiplier on grid size. nth-grid will be of size
+  # (initial_grid^grid_progression)^dim
+  set grid_progression  = 1.5
+
+  # Initial grid of size (initial_grid_size)^dim
+  set initial_grid_size = 3
+
+  # Number of grids in grid study
+  set number_of_grids   = 6
+end
diff --git a/tests/euler_implicit/2d_euler_implicit.prm b/tests/euler_implicit/2d_euler_implicit.prm
@@ -0,0 +1,41 @@
+# Listing of Parameters
+# ---------------------
+# Number of dimensions
+set dimension = 2
+
+# The PDE we want to solve. Choices are
+# <advection|diffusion|convection_diffusion>.
+set pde_type  = euler
+
+subsection ODE solver
+  # Maximum nonlinear solver iterations
+  set nonlinear_max_iterations            = 500000
+
+  # Nonlinear solver residual tolerance
+  set nonlinear_steady_residual_tolerance = 1e-13
+
+  # Print every print_iteration_modulo iterations of the nonlinear solver
+  set print_iteration_modulo              = 1
+
+  # Explicit or implicit solverChoices are <explicit|implicit>.
+  set ode_solver_type                         = implicit
+end
+
+subsection manufactured solution convergence study
+  # Last degree used for convergence study
+  set degree_end        = 3
+
+  # Starting degree for convergence study
+  set degree_start      = 1
+
+  # Multiplier on grid size. nth-grid will be of size
+  # (initial_grid^grid_progression)^dim
+  set grid_progression  = 1.5
+
+  # Initial grid of size (initial_grid_size)^dim
+  set initial_grid_size = 2
+
+  # Number of grids in grid study
+  set number_of_grids   = 5
+end
+
diff --git a/tests/euler_implicit/3d_euler_implicit.prm b/tests/euler_implicit/3d_euler_implicit.prm
@@ -0,0 +1,42 @@
+# Listing of Parameters
+# ---------------------
+# Number of dimensions
+set dimension = 3
+
+# The PDE we want to solve. Choices are
+# <advection|diffusion|convection_diffusion>.
+set pde_type  = euler
+
+subsection ODE solver
+  # Maximum nonlinear solver iterations
+  set nonlinear_max_iterations            = 500000
+
+  # Nonlinear solver residual tolerance
+  set nonlinear_steady_residual_tolerance = 1e-13
+
+  # Print every print_iteration_modulo iterations of the nonlinear solver
+  set print_iteration_modulo              = 1
+
+  # Explicit or implicit solverChoices are <explicit|implicit>.
+  set ode_solver_type                         = implicit
+end
+
+subsection manufactured solution convergence study
+  # Last degree used for convergence study
+  set degree_end        = 3
+
+  # Starting degree for convergence study
+  set degree_start      = 1
+
+  # Multiplier on grid size. nth-grid will be of size
+  # (initial_grid^grid_progression)^dim
+  set grid_progression  = 1.5
+
+  # Initial grid of size (initial_grid_size)^dim
+  set initial_grid_size = 2
+
+  # Number of grids in grid study
+  set number_of_grids   = 4
+end
+
+
diff --git a/tests/euler_implicit/CMakeLists.txt b/tests/euler_implicit/CMakeLists.txt
@@ -0,0 +1,22 @@
+set(TEST_OUTPUT_DIR ${CMAKE_CURRENT_BINARY_DIR})
+
+configure_file(1d_euler_implicit.prm 1d_euler_implicit.prm COPYONLY)
+add_test(
+  NAME 1D_EULER_IMPLICIT_MANUFACTURED_SOLUTION
+  COMMAND ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_1D -i ${CMAKE_CURRENT_BINARY_DIR}/1d_euler_implicit.prm
+  WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
+)
+
+configure_file(2d_euler_implicit.prm 2d_euler_implicit.prm COPYONLY)
+add_test(
+  NAME 2D_EULER_IMPLICIT_MANUFACTURED_SOLUTION
+  COMMAND ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_2D -i ${CMAKE_CURRENT_BINARY_DIR}/2d_euler_implicit.prm
+  WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
+)
+
+configure_file(3d_euler_implicit.prm 3d_euler_implicit.prm COPYONLY)
+add_test(
+  NAME 3D_EULER_IMPLICIT_MANUFACTURED_SOLUTION
+  COMMAND ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_3D -i ${CMAKE_CURRENT_BINARY_DIR}/3d_euler_implicit.prm
+  WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
+)