From c256b26a0700b9854d9cb6ad25ea012b8ae964ff Mon Sep 17 00:00:00 2001 From: "A. Cicchino" <39313421+AlexanderCicchino@users.noreply.github.com> Date: Sat, 27 Jan 2024 18:34:22 -0500 Subject: [PATCH] Fixed exterior reference normal zone issue in strong. (#250) This PR fixes/closes the zone issue in strong DG as per issue #249 . The problem was that previously, on a face, it was assumed that the reference faces were neighbours, which is the case for grids made from hypercube and a warping or a single mesh in dealii. In the general case, this isn't, and now it has been fixed. Note, the reference cell itself could be rotated and no longer be a square going from 0->1, that case has not been tested in strong. ![strong_naca0012_solution](https://github.com/dougshidong/PHiLiP/assets/39313421/102c466b-ecc6-4268-90d0-17a7557650ff) --------- Co-authored-by: Carolyn --- src/dg/strong_dg.cpp | 40 ++++--- src/flow_solver/flow_solver_cases/naca0012.h | 2 +- src/parameters/all_parameters.cpp | 3 + src/parameters/all_parameters.h | 1 + src/testing/CMakeLists.txt | 1 + src/testing/naca0012_unsteady_check_quick.cpp | 113 ++++++++++++++++++ src/testing/naca0012_unsteady_check_quick.h | 40 +++++++ src/testing/tests.cpp | 3 + .../euler_integration/naca0012/CMakeLists.txt | 7 ++ .../naca0012/naca0012_unsteady.prm | 36 ++++++ 10 files changed, 226 insertions(+), 20 deletions(-) create mode 100644 src/testing/naca0012_unsteady_check_quick.cpp create mode 100644 src/testing/naca0012_unsteady_check_quick.h create mode 100644 tests/integration_tests_control_files/euler_integration/naca0012/naca0012_unsteady.prm diff --git a/src/dg/strong_dg.cpp b/src/dg/strong_dg.cpp index 06dc710b0..b08255c6e 100644 --- a/src/dg/strong_dg.cpp +++ b/src/dg/strong_dg.cpp @@ -2096,8 +2096,10 @@ void DGStrong::assemble_face_term_strong( // we exploit the fact that the unit reference normal has a value of 0 in all reference directions except // the outward reference normal dircetion. const dealii::Tensor<1,dim,double> unit_ref_normal_int = dealii::GeometryInfo::unit_normal_vector[iface]; + const dealii::Tensor<1,dim,double> unit_ref_normal_ext = dealii::GeometryInfo::unit_normal_vector[neighbor_iface]; // Extract the reference direction that is outward facing on the facet. - const int dim_not_zero = iface / 2;//reference direction of face integer division + const int dim_not_zero_int = iface / 2;//reference direction of face integer division + const int dim_not_zero_ext = neighbor_iface / 2;//reference direction of face integer division std::array,nstate> conv_int_vol_ref_flux_interp_to_face_dot_ref_normal; std::array,nstate> conv_ext_vol_ref_flux_interp_to_face_dot_ref_normal; @@ -2116,32 +2118,32 @@ void DGStrong::assemble_face_term_strong( // interpolate reference volume convective flux to the facet, and apply unit reference normal as scaled by 1.0 or -1.0 if(!this->all_parameters->use_split_form && !this->all_parameters->use_curvilinear_split_form){ flux_basis_int.matrix_vector_mult_surface_1D(iface, - conv_ref_flux_at_vol_q_int[istate][dim_not_zero], + conv_ref_flux_at_vol_q_int[istate][dim_not_zero_int], conv_int_vol_ref_flux_interp_to_face_dot_ref_normal[istate], flux_basis_int.oneD_surf_operator,//the flux basis interpolates from the flux nodes flux_basis_int.oneD_vol_operator, - false, unit_ref_normal_int[dim_not_zero]);//don't add to previous value, scale by unit_normal int + false, unit_ref_normal_int[dim_not_zero_int]);//don't add to previous value, scale by unit_normal int flux_basis_ext.matrix_vector_mult_surface_1D(neighbor_iface, - conv_ref_flux_at_vol_q_ext[istate][dim_not_zero], + conv_ref_flux_at_vol_q_ext[istate][dim_not_zero_ext], conv_ext_vol_ref_flux_interp_to_face_dot_ref_normal[istate], flux_basis_ext.oneD_surf_operator, flux_basis_ext.oneD_vol_operator, - false, -unit_ref_normal_int[dim_not_zero]);//don't add to previous value, unit_normal ext is -unit normal int + false, unit_ref_normal_ext[dim_not_zero_ext]);//don't add to previous value, unit_normal ext is -unit normal int } // interpolate reference volume dissipative flux to the facet, and apply unit reference normal as scaled by 1.0 or -1.0 flux_basis_int.matrix_vector_mult_surface_1D(iface, - diffusive_ref_flux_at_vol_q_int[istate][dim_not_zero], + diffusive_ref_flux_at_vol_q_int[istate][dim_not_zero_int], diffusive_int_vol_ref_flux_interp_to_face_dot_ref_normal[istate], flux_basis_int.oneD_surf_operator, flux_basis_int.oneD_vol_operator, - false, unit_ref_normal_int[dim_not_zero]); + false, unit_ref_normal_int[dim_not_zero_int]); flux_basis_ext.matrix_vector_mult_surface_1D(neighbor_iface, - diffusive_ref_flux_at_vol_q_ext[istate][dim_not_zero], + diffusive_ref_flux_at_vol_q_ext[istate][dim_not_zero_ext], diffusive_ext_vol_ref_flux_interp_to_face_dot_ref_normal[istate], flux_basis_ext.oneD_surf_operator, flux_basis_ext.oneD_vol_operator, - false, -unit_ref_normal_int[dim_not_zero]); + false, unit_ref_normal_ext[dim_not_zero_ext]); } @@ -2235,8 +2237,8 @@ void DGStrong::assemble_face_term_strong( std::vector Hadamard_rows_sparsity_ext(row_size_ext); std::vector Hadamard_columns_sparsity_ext(col_size_ext); if(this->all_parameters->use_split_form || this->all_parameters->use_curvilinear_split_form){ - flux_basis_int.sum_factorized_Hadamard_surface_sparsity_pattern(n_face_quad_pts, n_quad_pts_1D_int, Hadamard_rows_sparsity_int, Hadamard_columns_sparsity_int, dim_not_zero); - flux_basis_ext.sum_factorized_Hadamard_surface_sparsity_pattern(n_face_quad_pts, n_quad_pts_1D_ext, Hadamard_rows_sparsity_ext, Hadamard_columns_sparsity_ext, dim_not_zero); + flux_basis_int.sum_factorized_Hadamard_surface_sparsity_pattern(n_face_quad_pts, n_quad_pts_1D_int, Hadamard_rows_sparsity_int, Hadamard_columns_sparsity_int, dim_not_zero_int); + flux_basis_ext.sum_factorized_Hadamard_surface_sparsity_pattern(n_face_quad_pts, n_quad_pts_1D_ext, Hadamard_rows_sparsity_ext, Hadamard_columns_sparsity_ext, dim_not_zero_ext); } std::array,nstate> surf_vol_ref_2pt_flux_interp_surf_int; @@ -2313,7 +2315,7 @@ void DGStrong::assemble_face_term_strong( 0.5*(metric_cofactor_surf + metric_cofactor_vol_int), conv_ref_flux_2pt); //only store the dim not zero in reference space bc dot product with unit ref normal later. - surface_ref_2pt_flux_int[istate][iquad_face][column_index] = conv_ref_flux_2pt[dim_not_zero]; + surface_ref_2pt_flux_int[istate][iquad_face][column_index] = conv_ref_flux_2pt[dim_not_zero_int]; } } for(unsigned int row_index = iquad_face * n_quad_pts_1D_ext, column_index = 0; @@ -2352,7 +2354,7 @@ void DGStrong::assemble_face_term_strong( 0.5*(metric_cofactor_surf + metric_cofactor_vol_ext), conv_ref_flux_2pt); //only store the dim not zero in reference space bc dot product with unit ref normal later. - surface_ref_2pt_flux_ext[istate][iquad_face][column_index] = conv_ref_flux_2pt[dim_not_zero]; + surface_ref_2pt_flux_ext[istate][iquad_face][column_index] = conv_ref_flux_2pt[dim_not_zero_ext]; } } } @@ -2368,7 +2370,7 @@ void DGStrong::assemble_face_term_strong( flux_basis_int.oneD_surf_operator[iface_1D], oneD_quad_weights_vol_int, surf_oper_sparse_int, - dim_not_zero); + dim_not_zero_int); const int neighbor_iface_1D = neighbor_iface % 2;//the reference neighbour face number const std::vector &oneD_quad_weights_vol_ext = this->oneD_quadrature_collection[poly_degree_ext].get_weights(); dealii::FullMatrix surf_oper_sparse_ext(n_face_quad_pts, n_quad_pts_1D_ext); @@ -2377,7 +2379,7 @@ void DGStrong::assemble_face_term_strong( flux_basis_ext.oneD_surf_operator[neighbor_iface_1D], oneD_quad_weights_vol_ext, surf_oper_sparse_ext, - dim_not_zero); + dim_not_zero_ext); // Apply the surface Hadamard products and multiply with vector of ones for both off diagonal terms in // Eq.(15) in Chan, Jesse. "Skew-symmetric entropy stable modal discontinuous Galerkin formulations." Journal of Scientific Computing 81.1 (2019): 459-485. @@ -2401,20 +2403,20 @@ void DGStrong::assemble_face_term_strong( for(unsigned int iquad_int=0; iquad_int const std::shared_ptr > dg, const std::shared_ptr unsteady_data_table) override; -private: +public: /// Compute lift double compute_lift(std::shared_ptr> dg) const; diff --git a/src/parameters/all_parameters.cpp b/src/parameters/all_parameters.cpp index ff1925a10..b3d34b22a 100644 --- a/src/parameters/all_parameters.cpp +++ b/src/parameters/all_parameters.cpp @@ -184,6 +184,7 @@ void AllParameters::declare_parameters (dealii::ParameterHandler &prm) " burgers_energy_conservation_rrk | " " euler_entropy_conserving_split_forms_check | " " h_refinement_study_isentropic_vortex | " + " naca0012_unsteady_check_quick | " " khi_robustness"), "The type of test we want to solve. " "Choices are " @@ -222,6 +223,7 @@ void AllParameters::declare_parameters (dealii::ParameterHandler &prm) " burgers_energy_conservation_rrk | " " euler_entropy_conserving_split_forms_check | " " h_refinement_study_isentropic_vortex | " + " naca0012_unsteady_check_quick | " " khi_robustness>."); prm.declare_entry("pde_type", "advection", @@ -398,6 +400,7 @@ void AllParameters::parse_parameters (dealii::ParameterHandler &prm) { test_type = euler_entropy_conserving_split_forms_check; } else if (test_string == "h_refinement_study_isentropic_vortex") { test_type = h_refinement_study_isentropic_vortex; } else if (test_string == "khi_robustness") { test_type = khi_robustness; } + else if (test_string == "naca0012_unsteady_check_quick") { test_type = naca0012_unsteady_check_quick; } overintegration = prm.get_integer("overintegration"); diff --git a/src/parameters/all_parameters.h b/src/parameters/all_parameters.h index 70cd270a9..9a07bb46e 100644 --- a/src/parameters/all_parameters.h +++ b/src/parameters/all_parameters.h @@ -192,6 +192,7 @@ class AllParameters euler_entropy_conserving_split_forms_check, h_refinement_study_isentropic_vortex, khi_robustness, + naca0012_unsteady_check_quick, homogeneous_isotropic_turbulence_initialization_check, }; /// Store selected TestType from the input file. diff --git a/src/testing/CMakeLists.txt b/src/testing/CMakeLists.txt index 6d13c791d..7de9c30d2 100644 --- a/src/testing/CMakeLists.txt +++ b/src/testing/CMakeLists.txt @@ -35,6 +35,7 @@ set(TESTS_SOURCE euler_entropy_conserving_split_forms_check.cpp homogeneous_isotropic_turbulence_initialization_check.cpp khi_robustness.cpp + naca0012_unsteady_check_quick.cpp ) foreach(dim RANGE 1 3) diff --git a/src/testing/naca0012_unsteady_check_quick.cpp b/src/testing/naca0012_unsteady_check_quick.cpp new file mode 100644 index 000000000..45ec426e6 --- /dev/null +++ b/src/testing/naca0012_unsteady_check_quick.cpp @@ -0,0 +1,113 @@ +#include "naca0012_unsteady_check_quick.h" +#include "flow_solver/flow_solver_factory.h" +#include "flow_solver/flow_solver_cases/naca0012.h" + +namespace PHiLiP { +namespace Tests { + +template +NACA0012UnsteadyCheckQuick::NACA0012UnsteadyCheckQuick( + const PHiLiP::Parameters::AllParameters *const parameters_input, + const dealii::ParameterHandler ¶meter_handler_input) + : TestsBase::TestsBase(parameters_input) + , parameter_handler(parameter_handler_input) +{} + +template +Parameters::AllParameters NACA0012UnsteadyCheckQuick::reinit_params(const bool use_weak_form_input, const bool use_two_point_flux_input) const +{ + PHiLiP::Parameters::AllParameters parameters = *(this->all_parameters); + + parameters.use_weak_form = use_weak_form_input; + using ConvFluxEnum = Parameters::AllParameters::ConvectiveNumericalFlux; + parameters.conv_num_flux_type = (use_two_point_flux_input) ? ConvFluxEnum::two_point_flux : ConvFluxEnum::roe; + + return parameters; +} + +template +int NACA0012UnsteadyCheckQuick::run_test() const +{ + const int n_runs = 3; + double lift_calculated[n_runs] = {0}; + double drag_calculated[n_runs] = {0}; + + const bool use_weak_form[3] = {true, false, false}; // choose weak or strong + const bool use_two_point_flux[3] = {false, false, true}; // choose two point flux or roe flux + for (unsigned int irun = 0; irun < n_runs; ++irun) { + + this->pcout << "=====================================================" << std::endl; + // Make new parameters according to the current run + const Parameters::AllParameters params_loop = reinit_params(use_weak_form[irun], use_two_point_flux[irun]); + + if (use_weak_form[irun]){ + this->pcout << "Initialized parameters with weak form." << std::endl; + } else{ + this->pcout << "Initialized parameters with strong form." << std::endl; + } + + if (use_two_point_flux[irun]){ + this->pcout << "Using two-point flux." << std::endl; + } else{ + this->pcout << "Using roe flux." << std::endl; + } + this->pcout << std::endl; + + // Initialize flow_solver + std::unique_ptr> flow_solver_loop = FlowSolver::FlowSolverFactory::select_flow_case(¶ms_loop, parameter_handler); + + static_cast(flow_solver_loop->run()); + std::unique_ptr> flow_solver_case_loop = std::make_unique>(this->all_parameters); + lift_calculated[irun] = flow_solver_case_loop->compute_lift((flow_solver_loop->dg)); + drag_calculated[irun] = flow_solver_case_loop->compute_drag((flow_solver_loop->dg)); + this->pcout << "Finished run." << std::endl; + this->pcout << "Calculated lift value was : " << lift_calculated[irun] << std::endl + << "Calculated drag value was : " << drag_calculated[irun] << std::endl; + this->pcout << "=====================================================" << std::endl; + } + + const double acceptable_tolerance = 0.00001; + int testfail = 0; + + // For now, allow split form to have a different end value. Leaving as a condition so we can reevaluate this choice later. + const bool allow_strong_split_to_be_different = true; + + this->pcout << std::endl + << "Finished runs. Summary of results: " << std::endl + << "Scheme | Lift | Drag" << std::endl + << "------------------------------------" << std::endl + << "Weak, roe | " << lift_calculated[0] << " | " << drag_calculated[0] << std::endl + << "Strong, roe | " << lift_calculated[1] << " | " << drag_calculated[1] << std::endl + << "Strong, split | " << lift_calculated[2] << " | " << drag_calculated[2] << std::endl; + + if (allow_strong_split_to_be_different) { + if ((abs(lift_calculated[0]-lift_calculated[1]) > acceptable_tolerance) + || (abs(drag_calculated[0]-drag_calculated[1]) > acceptable_tolerance)){ + testfail = 1; + } + } else{ + const double lift_max = *std::max_element(lift_calculated, lift_calculated+n_runs); + const double lift_min = *std::min_element(lift_calculated, lift_calculated+n_runs); + const double drag_max = *std::max_element(drag_calculated, drag_calculated+n_runs); + const double drag_min = *std::min_element(drag_calculated, drag_calculated+n_runs); + if ((abs(lift_max - lift_min) > acceptable_tolerance) + || (abs(drag_max - drag_min) > acceptable_tolerance)){ + testfail = 1; + } + + } + if (testfail == 1){ + this->pcout << "Test failing: difference is not within the allowable tolerance of " << acceptable_tolerance << std::endl; + this->pcout << "If this test is failing, please check the *vtu output." << std::endl; + } else { + this->pcout << "Test passing: lift and drag close to the expected value." << std::endl; + } + + return testfail; +} + +#if PHILIP_DIM==2 + template class NACA0012UnsteadyCheckQuick; +#endif +} // Tests namespace +} // PHiLiP namespace diff --git a/src/testing/naca0012_unsteady_check_quick.h b/src/testing/naca0012_unsteady_check_quick.h new file mode 100644 index 000000000..fe8295ea6 --- /dev/null +++ b/src/testing/naca0012_unsteady_check_quick.h @@ -0,0 +1,40 @@ + +#ifndef __NACA0012_UNSTEADY_CHECK_QUICK__ +#define __NACA0012_UNSTEADY_CHECK_QUICK__ + +#include "tests.h" + +namespace PHiLiP { +namespace Tests { + +/// NACA 0012 Unsteady Check +/** Runs a short time interval for unsteady Euler flow over NACA0012 airfoil + * Uses GMSH reader for the mesh + * Compares against a hard-coded expectation value for lift + * NOTE: it has not been verified that the results are physically meaningful; + * this is a verification test that weak and strong give relatively consistent results. + */ +template +class NACA0012UnsteadyCheckQuick: public TestsBase +{ +public: + /// Constructor + NACA0012UnsteadyCheckQuick( + const Parameters::AllParameters *const parameters_input, + const dealii::ParameterHandler ¶meter_handler_input); + + /// Parameter handler for storing the .prm file being ran + const dealii::ParameterHandler ¶meter_handler; + + /// Run test + int run_test () const override; +protected: + + /// Reinit parameters based on specified weak/strong form and convective num flux + Parameters::AllParameters reinit_params(const bool use_weak_form_input, const bool use_two_point_flux_input) const; +}; + +} // End of Tests namespace +} // End of PHiLiP namespace + +#endif diff --git a/src/testing/tests.cpp b/src/testing/tests.cpp index 5c2898bfd..a3652ea9f 100644 --- a/src/testing/tests.cpp +++ b/src/testing/tests.cpp @@ -43,6 +43,7 @@ #include "euler_entropy_conserving_split_forms_check.h" #include "homogeneous_isotropic_turbulence_initialization_check.h" #include "khi_robustness.h" +#include "naca0012_unsteady_check_quick.h" namespace PHiLiP { namespace Tests { @@ -301,6 +302,8 @@ ::select_test(const AllParam *const parameters_input, if constexpr (dim==3 && nstate==dim+2) return std::make_unique>(parameters_input, parameter_handler_input); } else if(test_type == Test_enum::khi_robustness) { if constexpr (dim==2 && nstate==dim+2) return std::make_unique>(parameters_input, parameter_handler_input); + } else if(test_type == Test_enum::naca0012_unsteady_check_quick){ + if constexpr (dim==2 && nstate==dim+2) return std::make_unique>(parameters_input, parameter_handler_input); } else { std::cout << "Invalid test. You probably forgot to add it to the list of tests in tests.cpp" << std::endl; std::abort(); diff --git a/tests/integration_tests_control_files/euler_integration/naca0012/CMakeLists.txt b/tests/integration_tests_control_files/euler_integration/naca0012/CMakeLists.txt index c612c58cc..5339c0eea 100644 --- a/tests/integration_tests_control_files/euler_integration/naca0012/CMakeLists.txt +++ b/tests/integration_tests_control_files/euler_integration/naca0012/CMakeLists.txt @@ -33,6 +33,13 @@ add_test( WORKING_DIRECTORY ${TEST_OUTPUT_DIR} ) +configure_file(naca0012_unsteady.prm naca0012_unsteady.prm COPYONLY) +add_test( + NAME 2D_EULER_INTEGRATION_UNSTEADY_NACA0012_SUBSONIC_STRONGvWEAK + COMMAND mpirun -np ${MPIMAX} ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_2D -i ${CMAKE_CURRENT_BINARY_DIR}/naca0012_unsteady.prm + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} +) + #configure_file(2d_navier_stokes_naca0012_subsonic_05_200.prm 2d_navier_stokes_naca0012_subsonic_05_200.prm COPYONLY) #add_test( # NAME 2D_NAVIER_STOKES_INTEGRATION_NACA0012_SUBSONIC_LONG diff --git a/tests/integration_tests_control_files/euler_integration/naca0012/naca0012_unsteady.prm b/tests/integration_tests_control_files/euler_integration/naca0012/naca0012_unsteady.prm new file mode 100644 index 000000000..bb3f9d800 --- /dev/null +++ b/tests/integration_tests_control_files/euler_integration/naca0012/naca0012_unsteady.prm @@ -0,0 +1,36 @@ +# Listing of Parameters +# --------------------- + +set test_type = naca0012_unsteady_check_quick +set dimension = 2 + +set use_weak_form = true # will be changed by test [strong/weak] +set pde_type = euler +set conv_num_flux = roe # will be changed by test [roe/two_point_flux] +set overintegration = 2 +set two_point_num_flux_type = Ra # only used if conv_num_flux == two_point_flux + +subsection euler + set reference_length = 1.0 + set mach_infinity = 0.50 + set angle_of_attack = 2.0 +end + +subsection ODE solver + set output_solution_every_x_steps = 10 + set ode_solver_type = runge_kutta + set print_iteration_modulo = 1 +end + +subsection flow_solver + set flow_case_type = naca0012 + set steady_state = false + set poly_degree = 3 + set constant_time_step = 0.00005 # no testing has been done to ensure that this is stable for many time steps; it is only guaranteed to reach final_time as set here. + set final_time = 0.0005 + subsection grid + set input_mesh_filename = ../../../meshes/naca0012_hopw_ref2 + end +end + +