1 #include "positivity_preserving_limiter.h" 2 #include "tvb_limiter.h" 11 template <
int dim,
int nstate,
typename real>
18 PDE_enum pde_type = parameters_input->
pde_type;
20 std::shared_ptr< ManufacturedSolutionFunction<dim, real> > manufactured_solution_function
23 if (pde_type == PDE_enum::euler &&
nstate == dim + 2) {
24 euler_physics = std::make_shared < Physics::Euler<dim, nstate, real> >(
27 parameters_input->euler_param.gamma_gas,
28 parameters_input->euler_param.mach_inf,
29 parameters_input->euler_param.angle_of_attack,
30 parameters_input->euler_param.side_slip_angle,
31 manufactured_solution_function,
32 parameters_input->two_point_num_flux_type);
35 std::cout <<
"Error: Positivity-Preserving Limiter can only be applied for pde_type==euler" << std::endl;
42 tvbLimiter = std::make_shared < TVBLimiter<dim, nstate, real> >(parameters_input);
45 std::cout <<
"Error: Cannot create TVB limiter for dim > 1" << std::endl;
51 template <
int dim,
int nstate,
typename real>
53 const std::vector< real >& p_lim,
54 const std::array<real, nstate>& soln_cell_avg,
55 const std::array<std::vector<real>,
nstate>& soln_at_q,
56 const unsigned int n_quad_pts,
60 std::vector<real> theta2(n_quad_pts, 1);
61 for (
unsigned int iquad = 0; iquad < n_quad_pts; ++iquad) {
62 if (p_lim[iquad] >= eps) {
66 real s_coeff1 = (soln_at_q[
nstate - 1][iquad] - soln_cell_avg[
nstate - 1]) * (soln_at_q[0][iquad] - soln_cell_avg[0]) - 0.5 * pow(soln_at_q[1][iquad] - soln_cell_avg[1], 2);
68 real s_coeff2 = soln_cell_avg[0] * (soln_at_q[
nstate - 1][iquad] - soln_cell_avg[
nstate - 1]) + soln_cell_avg[
nstate - 1] * (soln_at_q[0][iquad] - soln_cell_avg[0])
69 - soln_cell_avg[1] * (soln_at_q[1][iquad] - soln_cell_avg[1]) - (eps / gamma) * (soln_at_q[0][iquad] - soln_cell_avg[0]);
71 real s_coeff3 = (soln_cell_avg[
nstate - 1] * soln_cell_avg[0]) - 0.5 * pow(soln_cell_avg[1], 2) - (eps / gamma) * soln_cell_avg[0];
74 s_coeff1 -= 0.5 * pow(soln_at_q[2][iquad] - soln_cell_avg[2], 2);
75 s_coeff2 -= soln_cell_avg[2] * (soln_at_q[2][iquad] - soln_cell_avg[2]);
76 s_coeff3 -= 0.5 * pow(soln_cell_avg[2], 2);
80 s_coeff1 -= 0.5 * pow(soln_at_q[3][iquad] - soln_cell_avg[3], 2);
81 s_coeff2 -= soln_cell_avg[3] * (soln_at_q[3][iquad] - soln_cell_avg[3]);
82 s_coeff3 -= 0.5 * pow(soln_cell_avg[3], 2);
85 real discriminant = s_coeff2 * s_coeff2 - 4 * s_coeff1 * s_coeff3;
87 if (discriminant >= 0) {
88 real t1 = (-s_coeff2 + sqrt(discriminant)) / (2 * s_coeff1);
89 real t2 = (-s_coeff2 - sqrt(discriminant)) / (2 * s_coeff1);
90 theta2[iquad] = std::min(t1, t2);
101 template <
int dim,
int nstate,
typename real>
103 const std::array<std::vector<real>,
nstate>& soln_at_q,
104 const unsigned int n_quad_pts,
107 std::vector<real> t2(n_quad_pts, 1);
109 std::array<real, nstate> soln_at_iquad;
112 for (
unsigned int iquad = 0; iquad < n_quad_pts; ++iquad) {
113 for (
unsigned int istate = 0; istate <
nstate; ++istate) {
114 soln_at_iquad[istate] = soln_at_q[istate][iquad];
118 if (nstate == dim + 2)
124 t2[iquad] = p_avg / (p_avg - p_lim);
126 if (t2[iquad] != 1) {
127 if (t2[iquad] >= 0 && t2[iquad] <= 1 && t2[iquad] < theta2) {
136 template <
int dim,
int nstate,
typename real>
138 const double density_avg,
139 const double density_min,
140 const double lower_bound,
144 real eps = std::min({ lower_bound, density_avg, p_avg });
145 if (eps < 0) eps = lower_bound;
148 if (density_avg - density_min > 1e-13)
149 theta = (density_avg - density_min == 0) ? 1.0 : std::min((density_avg - eps) / (density_avg - density_min), 1.0);
151 if (theta < 0 || theta > 1)
157 template <
int dim,
int nstate,
typename real>
159 dealii::LinearAlgebra::distributed::Vector<double>& solution,
160 const std::array<std::vector<real>,
nstate>& soln_coeff,
161 const unsigned int n_shape_fns,
162 const std::vector<dealii::types::global_dof_index>& current_dofs_indices)
165 for (
int istate = 0; istate <
nstate; istate++) {
166 for (
unsigned int ishape = 0; ishape < n_shape_fns; ++ishape) {
167 const unsigned int idof = istate * n_shape_fns + ishape;
168 solution[current_dofs_indices[idof]] = soln_coeff[istate][ishape];
171 if (istate == 0 && solution[current_dofs_indices[idof]] < 0) {
172 std::cout <<
"Error: Density is a negative value - Aborting... " << std::endl << solution[current_dofs_indices[idof]] << std::endl << std::flush;
179 template <
int dim,
int nstate,
typename real>
181 dealii::LinearAlgebra::distributed::Vector<double>& solution,
182 const dealii::DoFHandler<dim>& dof_handler,
183 const dealii::hp::FECollection<dim>& fe_collection,
184 const dealii::hp::QCollection<dim>& volume_quadrature_collection,
185 const unsigned int grid_degree,
186 const unsigned int max_degree,
187 const dealii::hp::FECollection<1> oneD_fe_collection_1state,
188 const dealii::hp::QCollection<1> oneD_quadrature_collection)
192 this->
tvbLimiter->limit(solution, dof_handler, fe_collection, volume_quadrature_collection, grid_degree, max_degree, oneD_fe_collection_1state, oneD_quadrature_collection);
197 const unsigned int init_grid_degree = grid_degree;
205 soln_basis_projection_oper.
build_1D_volume_operator(oneD_fe_collection_1state[max_degree], oneD_quadrature_collection[max_degree]);
207 for (
auto soln_cell : dof_handler.active_cell_iterators()) {
208 if (!soln_cell->is_locally_owned())
continue;
210 std::vector<dealii::types::global_dof_index> current_dofs_indices;
212 const int i_fele = soln_cell->active_fe_index();
213 const dealii::FESystem<dim, dim>& current_fe_ref = fe_collection[i_fele];
214 const int poly_degree = current_fe_ref.tensor_degree();
216 const unsigned int n_dofs_curr_cell = current_fe_ref.n_dofs_per_cell();
219 current_dofs_indices.resize(n_dofs_curr_cell);
220 soln_cell->get_dof_indices(current_dofs_indices);
223 std::array<std::vector<real>,
nstate> soln_coeff;
225 const unsigned int n_shape_fns = n_dofs_curr_cell /
nstate;
226 real local_min_density = 1e6;
228 for (
unsigned int istate = 0; istate <
nstate; ++istate) {
229 soln_coeff[istate].resize(n_shape_fns);
233 for (
unsigned int idof = 0; idof < n_dofs_curr_cell; ++idof) {
234 const unsigned int istate = fe_collection[poly_degree].system_to_component_index(idof).first;
235 const unsigned int ishape = fe_collection[poly_degree].system_to_component_index(idof).second;
236 soln_coeff[istate][ishape] = solution[current_dofs_indices[idof]];
242 const unsigned int n_quad_pts = volume_quadrature_collection[poly_degree].size();
243 const std::vector<real>& quad_weights = volume_quadrature_collection[poly_degree].get_weights();
244 std::array<std::vector<real>, nstate> soln_at_q;
247 for (
int istate = 0; istate <
nstate; istate++) {
248 soln_at_q[istate].resize(n_quad_pts);
253 for (
unsigned int iquad = 0; iquad < n_quad_pts; ++iquad) {
254 if (soln_at_q[0][iquad] < local_min_density)
255 local_min_density = soln_at_q[0][iquad];
258 std::array<real, nstate> soln_cell_avg =
get_soln_cell_avg(soln_at_q, n_quad_pts, quad_weights);
263 if (nstate == dim + 2) {
272 for (
unsigned int iquad = 0; iquad < n_quad_pts; ++iquad) {
273 if (isnan(soln_at_q[0][iquad])) {
274 std::cout <<
"Error: Density at quadrature point is NaN - Aborting... " << std::endl << std::flush;
277 soln_at_q[0][iquad] = theta * (soln_at_q[0][iquad] - soln_cell_avg[0])
284 if (limiter_type == limiter_enum::positivity_preservingZhang2010 && nstate == dim + 2) {
285 std::vector< real > p_lim(n_quad_pts);
286 std::array<real, nstate> soln_at_iquad;
289 for (
unsigned int iquad = 0; iquad < n_quad_pts; ++iquad) {
290 for (
unsigned int istate = 0; istate <
nstate; ++istate) {
291 soln_at_iquad[istate] = soln_at_q[istate][iquad];
293 p_lim[iquad] =
euler_physics->compute_pressure(soln_at_iquad);
300 for (
unsigned int istate = 0; istate <
nstate; ++istate) {
301 for (
unsigned int iquad = 0; iquad < n_quad_pts; ++iquad) {
302 soln_at_q[istate][iquad] = theta2[iquad] * (soln_at_q[istate][iquad] - soln_cell_avg[istate])
303 + soln_cell_avg[istate];
308 else if (limiter_type == limiter_enum::positivity_preservingWang2012 && nstate == dim + 2) {
312 for (
unsigned int istate = 0; istate <
nstate; ++istate) {
313 for (
unsigned int iquad = 0; iquad < n_quad_pts; ++iquad) {
314 soln_at_q[istate][iquad] = theta2 * (soln_at_q[istate][iquad] - soln_cell_avg[istate])
315 + soln_cell_avg[istate];
321 for (
int istate = 0; istate <
nstate; istate++) {
PartialDifferentialEquation pde_type
Store the PDE type to be solved.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
LimiterType
Limiter type to be applied on the solution.
LimiterParam limiter_param
Contains parameters for limiter.
real get_theta2_Wang2012(const std::array< std::vector< real >, nstate > &soln_at_q, const unsigned int n_quad_pts, const double p_avg)
Obtain the theta value used to scale all the states using 3.7 in Wang, Shu 2012.
LimiterType bound_preserving_limiter
Variable to store specified limiter type.
void write_limited_solution(dealii::LinearAlgebra::distributed::Vector< double > &solution, const std::array< std::vector< real >, nstate > &soln_coeff, const unsigned int n_shape_fns, const std::vector< dealii::types::global_dof_index > ¤t_dofs_indices)
std::shared_ptr< BoundPreservingLimiterState< dim, nstate, real > > tvbLimiter
Pointer to TVB limiter class (TVB limiter can be applied in conjunction with this limiter) ...
Base Class for bound preserving limiters templated on state.
double min_density
Epsilon value for Positivity-Preserving Limiter.
PositivityPreservingLimiter(const Parameters::AllParameters *const parameters_input)
Constructor.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
static std::shared_ptr< ManufacturedSolutionFunction< dim, real > > create_ManufacturedSolution(Parameters::AllParameters const *const param, int nstate)
Construct Manufactured solution object from global parameter file.
PartialDifferentialEquation
Possible Partial Differential Equations to solve.
const int nstate
Number of states.
Files for the baseline physics.
EulerParam euler_param
Contains parameters for the Euler equations non-dimensionalization.
Main parameter class that contains the various other sub-parameter classes.
std::array< real, nstate > get_soln_cell_avg(const std::array< std::vector< real >, nstate > &soln_at_q, const unsigned int n_quad_pts, const std::vector< real > &quad_weights)
Function to obtain the solution cell average.
std::vector< real > get_theta2_Zhang2010(const std::vector< real > &p_lim, const std::array< real, nstate > &soln_cell_avg, const std::array< std::vector< real >, nstate > &soln_at_q, const unsigned int n_quad_pts, const double eps, const double gamma)
Obtain the theta value used to scale all the states using 3.16-3.18 in Zhang, Shu 2010...
double ref_length
Reference length.
real get_density_scaling_value(const double density_avg, const double density_min, const double pos_eps, const double p_avg)
dealii::FullMatrix< double > oneD_vol_operator
Stores the one dimensional volume operator.
Class for implementation of two forms of the Positivity-Preserving limiter derived from BoundPreservi...
void limit(dealii::LinearAlgebra::distributed::Vector< double > &solution, const dealii::DoFHandler< dim > &dof_handler, const dealii::hp::FECollection< dim > &fe_collection, const dealii::hp::QCollection< dim > &volume_quadrature_collection, const unsigned int grid_degree, const unsigned int max_degree, const dealii::hp::FECollection< 1 > oneD_fe_collection_1state, const dealii::hp::QCollection< 1 > oneD_quadrature_collection)
bool use_tvb_limiter
Flag for applying TVB Limiter.
std::shared_ptr< Physics::Euler< dim, nstate, double > > euler_physics
Euler physics pointer. Used to compute pressure.
void matrix_vector_mult_1D(const std::vector< real > &input_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis_x, const bool adding=false, const double factor=1.0)
Apply the matrix vector operation using the 1D operator in each direction.
const Parameters::AllParameters *const all_parameters
Pointer to parameters object.
Projection operator corresponding to basis functions onto M-norm (L2).