Skip to content

Commit

Permalink
Fixed exterior reference normal zone issue in strong. (#250)
Browse files Browse the repository at this point in the history
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 <c.pethrick@gmail.com>
  • Loading branch information
AlexanderCicchino and cpethrick authored Jan 27, 2024
1 parent a29ecce commit c256b26
Show file tree
Hide file tree
Showing 10 changed files with 226 additions and 20 deletions.
40 changes: 21 additions & 19 deletions src/dg/strong_dg.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2096,8 +2096,10 @@ void DGStrong<dim,nstate,real,MeshType>::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<dim>::unit_normal_vector[iface];
const dealii::Tensor<1,dim,double> unit_ref_normal_ext = dealii::GeometryInfo<dim>::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<std::vector<real>,nstate> conv_int_vol_ref_flux_interp_to_face_dot_ref_normal;
std::array<std::vector<real>,nstate> conv_ext_vol_ref_flux_interp_to_face_dot_ref_normal;
Expand All @@ -2116,32 +2118,32 @@ void DGStrong<dim,nstate,real,MeshType>::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]);
}


Expand Down Expand Up @@ -2235,8 +2237,8 @@ void DGStrong<dim,nstate,real,MeshType>::assemble_face_term_strong(
std::vector<unsigned int> Hadamard_rows_sparsity_ext(row_size_ext);
std::vector<unsigned int> 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<std::vector<real>,nstate> surf_vol_ref_2pt_flux_interp_surf_int;
Expand Down Expand Up @@ -2313,7 +2315,7 @@ void DGStrong<dim,nstate,real,MeshType>::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;
Expand Down Expand Up @@ -2352,7 +2354,7 @@ void DGStrong<dim,nstate,real,MeshType>::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];
}
}
}
Expand All @@ -2368,7 +2370,7 @@ void DGStrong<dim,nstate,real,MeshType>::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<double> &oneD_quad_weights_vol_ext = this->oneD_quadrature_collection[poly_degree_ext].get_weights();
dealii::FullMatrix<real> surf_oper_sparse_ext(n_face_quad_pts, n_quad_pts_1D_ext);
Expand All @@ -2377,7 +2379,7 @@ void DGStrong<dim,nstate,real,MeshType>::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.
Expand All @@ -2401,20 +2403,20 @@ void DGStrong<dim,nstate,real,MeshType>::assemble_face_term_strong(
for(unsigned int iquad_int=0; iquad_int<n_quad_pts_1D_int; iquad_int++){
surf_vol_ref_2pt_flux_interp_surf_int[istate][iface_quad]
-= surface_ref_2pt_flux_int_Hadamard_with_surf_oper[iface_quad][iquad_int]
* unit_ref_normal_int[dim_not_zero];
* unit_ref_normal_int[dim_not_zero_int];
const unsigned int column_index = iface_quad * n_quad_pts_1D_int + iquad_int;
surf_vol_ref_2pt_flux_interp_vol_int[istate][Hadamard_columns_sparsity_int[column_index]]
+= surface_ref_2pt_flux_int_Hadamard_with_surf_oper[iface_quad][iquad_int]
* unit_ref_normal_int[dim_not_zero];
* unit_ref_normal_int[dim_not_zero_int];
}
for(unsigned int iquad_ext=0; iquad_ext<n_quad_pts_1D_ext; iquad_ext++){
surf_vol_ref_2pt_flux_interp_surf_ext[istate][iface_quad]
-= surface_ref_2pt_flux_ext_Hadamard_with_surf_oper[iface_quad][iquad_ext]
* (-unit_ref_normal_int[dim_not_zero]);
* (unit_ref_normal_ext[dim_not_zero_ext]);
const unsigned int column_index = iface_quad * n_quad_pts_1D_ext + iquad_ext;
surf_vol_ref_2pt_flux_interp_vol_ext[istate][Hadamard_columns_sparsity_ext[column_index]]
+= surface_ref_2pt_flux_ext_Hadamard_with_surf_oper[iface_quad][iquad_ext]
* (-unit_ref_normal_int[dim_not_zero]);
* (unit_ref_normal_ext[dim_not_zero_ext]);
}
}
}
Expand Down
2 changes: 1 addition & 1 deletion src/flow_solver/flow_solver_cases/naca0012.h
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ class NACA0012 : public FlowSolverCaseBase<dim,nstate>
const std::shared_ptr <DGBase<dim, double>> dg,
const std::shared_ptr<dealii::TableHandler> unsteady_data_table) override;

private:
public:
/// Compute lift
double compute_lift(std::shared_ptr<DGBase<dim, double>> dg) const;

Expand Down
3 changes: 3 additions & 0 deletions src/parameters/all_parameters.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 "
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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");

Expand Down
1 change: 1 addition & 0 deletions src/parameters/all_parameters.h
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
1 change: 1 addition & 0 deletions src/testing/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
113 changes: 113 additions & 0 deletions src/testing/naca0012_unsteady_check_quick.cpp
Original file line number Diff line number Diff line change
@@ -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 <int dim, int nstate>
NACA0012UnsteadyCheckQuick<dim, nstate>::NACA0012UnsteadyCheckQuick(
const PHiLiP::Parameters::AllParameters *const parameters_input,
const dealii::ParameterHandler &parameter_handler_input)
: TestsBase::TestsBase(parameters_input)
, parameter_handler(parameter_handler_input)
{}

template <int dim, int nstate>
Parameters::AllParameters NACA0012UnsteadyCheckQuick<dim,nstate>::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 dim, int nstate>
int NACA0012UnsteadyCheckQuick<dim, nstate>::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<FlowSolver::FlowSolver<dim,nstate>> flow_solver_loop = FlowSolver::FlowSolverFactory<dim,nstate>::select_flow_case(&params_loop, parameter_handler);

static_cast<void>(flow_solver_loop->run());
std::unique_ptr<FlowSolver::NACA0012<dim, nstate>> flow_solver_case_loop = std::make_unique<FlowSolver::NACA0012<dim,nstate>>(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<PHILIP_DIM,PHILIP_DIM+2>;
#endif
} // Tests namespace
} // PHiLiP namespace
40 changes: 40 additions & 0 deletions src/testing/naca0012_unsteady_check_quick.h
Original file line number Diff line number Diff line change
@@ -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 <int dim, int nstate>
class NACA0012UnsteadyCheckQuick: public TestsBase
{
public:
/// Constructor
NACA0012UnsteadyCheckQuick(
const Parameters::AllParameters *const parameters_input,
const dealii::ParameterHandler &parameter_handler_input);

/// Parameter handler for storing the .prm file being ran
const dealii::ParameterHandler &parameter_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
Loading

0 comments on commit c256b26

Please sign in to comment.