[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
positivity_preserving_limiter.cpp
1 #include "positivity_preserving_limiter.h"
2 #include "tvb_limiter.h"
3 
4 namespace PHiLiP {
5 /**********************************
6 *
7 * Positivity Preserving Limiter Class
8 *
9 **********************************/
10 // Constructor
11 template <int dim, int nstate, typename real>
13  const Parameters::AllParameters* const parameters_input)
14  : BoundPreservingLimiterState<dim,nstate,real>::BoundPreservingLimiterState(parameters_input)
15 {
16  // Create pointer to Euler Physics to compute pressure if pde_type==euler
18  PDE_enum pde_type = parameters_input->pde_type;
19 
20  std::shared_ptr< ManufacturedSolutionFunction<dim, real> > manufactured_solution_function
22 
23  if (pde_type == PDE_enum::euler && nstate == dim + 2) {
24  euler_physics = std::make_shared < Physics::Euler<dim, nstate, real> >(
25  parameters_input,
26  parameters_input->euler_param.ref_length,
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);
33  }
34  else {
35  std::cout << "Error: Positivity-Preserving Limiter can only be applied for pde_type==euler" << std::endl;
36  std::abort();
37  }
38 
39  // Create pointer to TVB Limiter class if use_tvb_limiter==true && dim == 1
40  if (parameters_input->limiter_param.use_tvb_limiter) {
41  if (dim == 1) {
42  tvbLimiter = std::make_shared < TVBLimiter<dim, nstate, real> >(parameters_input);
43  }
44  else {
45  std::cout << "Error: Cannot create TVB limiter for dim > 1" << std::endl;
46  std::abort();
47  }
48  }
49 }
50 
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,
57  const double eps,
58  const double gamma)
59 {
60  std::vector<real> theta2(n_quad_pts, 1); // Value used to linearly scale state variables
61  for (unsigned int iquad = 0; iquad < n_quad_pts; ++iquad) {
62  if (p_lim[iquad] >= eps) {
63  theta2[iquad] = 1;
64  }
65  else {
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);
67 
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]);
70 
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];
72 
73  if (dim > 1) {
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);
77  }
78 
79  if (dim > 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);
83  }
84 
85  real discriminant = s_coeff2 * s_coeff2 - 4 * s_coeff1 * s_coeff3;
86 
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);
91  }
92  else {
93  theta2[iquad] = 0;
94  }
95  }
96  }
97 
98  return theta2;
99 }
100 
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,
105  const double p_avg)
106 {
107  std::vector<real> t2(n_quad_pts, 1);
108  real theta2 = 1.0; // Value used to linearly scale state variables
109  std::array<real, nstate> soln_at_iquad;
110 
111  // Obtain theta2 value
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];
115  }
116  real p_lim = 0;
117 
118  if (nstate == dim + 2)
119  p_lim = euler_physics->compute_pressure(soln_at_iquad);
120 
121  if (p_lim >= 0)
122  t2[iquad] = 1;
123  else
124  t2[iquad] = p_avg / (p_avg - p_lim);
125 
126  if (t2[iquad] != 1) {
127  if (t2[iquad] >= 0 && t2[iquad] <= 1 && t2[iquad] < theta2) {
128  theta2 = t2[iquad];
129  }
130  }
131  }
132 
133  return theta2;
134 }
135 
136 template <int dim, int nstate, typename real>
138  const double density_avg,
139  const double density_min,
140  const double lower_bound,
141  const double p_avg)
142 {
143  // Get epsilon (lower bound for rho) for theta limiter
144  real eps = std::min({ lower_bound, density_avg, p_avg });
145  if (eps < 0) eps = lower_bound;
146 
147  real theta = 1.0; // Value used to linearly scale density
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);
150 
151  if (theta < 0 || theta > 1)
152  theta = 0;
153 
154  return theta;
155 }
156 
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)
163 {
164  // Write limited solution dofs to the global solution vector.
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]; //
169 
170  // Verify that positivity of density is preserved after application of theta2 limiter
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;
173  std::abort();
174  }
175  }
176  }
177 }
178 
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)
189 {
190  // If use_tvb_limiter is true, apply TVB limiter before applying positivity-preserving limiter
191  if (this->all_parameters->limiter_param.use_tvb_limiter == true)
192  this->tvbLimiter->limit(solution, dof_handler, fe_collection, volume_quadrature_collection, grid_degree, max_degree, oneD_fe_collection_1state, oneD_quadrature_collection);
193 
194  //create 1D solution polynomial basis functions and corresponding projection operator
195  //to interpolate the solution to the quadrature nodes, and to project it back to the
196  //modal coefficients.
197  const unsigned int init_grid_degree = grid_degree;
198 
199  // Constructor for the operators
200  OPERATOR::basis_functions<dim, 2 * dim, real> soln_basis(1, max_degree, init_grid_degree);
201  OPERATOR::vol_projection_operator<dim, 2 * dim, real> soln_basis_projection_oper(1, max_degree, init_grid_degree);
202 
203  // Build the oneD operator to perform interpolation/projection
204  soln_basis.build_1D_volume_operator(oneD_fe_collection_1state[max_degree], oneD_quadrature_collection[max_degree]);
205  soln_basis_projection_oper.build_1D_volume_operator(oneD_fe_collection_1state[max_degree], oneD_quadrature_collection[max_degree]);
206 
207  for (auto soln_cell : dof_handler.active_cell_iterators()) {
208  if (!soln_cell->is_locally_owned()) continue;
209 
210  std::vector<dealii::types::global_dof_index> current_dofs_indices;
211  // Current reference element related to this physical cell
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();
215 
216  const unsigned int n_dofs_curr_cell = current_fe_ref.n_dofs_per_cell();
217 
218  // Obtain the mapping from local dof indices to global dof indices
219  current_dofs_indices.resize(n_dofs_curr_cell);
220  soln_cell->get_dof_indices(current_dofs_indices);
221 
222  // Extract the local solution dofs in the cell from the global solution dofs
223  std::array<std::vector<real>, nstate> soln_coeff;
224 
225  const unsigned int n_shape_fns = n_dofs_curr_cell / nstate;
226  real local_min_density = 1e6;
227 
228  for (unsigned int istate = 0; istate < nstate; ++istate) {
229  soln_coeff[istate].resize(n_shape_fns);
230  }
231 
232  // Allocate solution dofs and set local min
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]]; //
237 
238  // if (istate == 0 && soln_coeff[istate][ishape] < local_min_density)
239  // local_min_density = soln_coeff[istate][ishape];
240  }
241 
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;
245 
246  // Interpolate solution dofs to quadrature pts.
247  for (int istate = 0; istate < nstate; istate++) {
248  soln_at_q[istate].resize(n_quad_pts);
249  soln_basis.matrix_vector_mult_1D(soln_coeff[istate], soln_at_q[istate],
250  soln_basis.oneD_vol_operator);
251  }
252 
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];
256  }
257  // Obtain solution cell average
258  std::array<real, nstate> soln_cell_avg = get_soln_cell_avg(soln_at_q, n_quad_pts, quad_weights);
259 
260  real lower_bound = this->all_parameters->limiter_param.min_density;
261  real p_avg = 1e-13;
262 
263  if (nstate == dim + 2) {
264  // Compute average value of pressure using soln_cell_avg
265  p_avg = euler_physics->compute_pressure(soln_cell_avg);
266  }
267 
268  // Obtain value used to linearly scale density
269  real theta = get_density_scaling_value(soln_cell_avg[0], local_min_density, lower_bound, p_avg);
270 
271  // Apply limiter on density values at quadrature points
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;
275  std::abort();
276  }
277  soln_at_q[0][iquad] = theta * (soln_at_q[0][iquad] - soln_cell_avg[0])
278  + soln_cell_avg[0];
279  }
280 
281  using limiter_enum = Parameters::LimiterParam::LimiterType;
282  limiter_enum limiter_type = this->all_parameters->limiter_param.bound_preserving_limiter;
283 
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;
287 
288  // Compute pressure at quadrature points
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];
292  }
293  p_lim[iquad] = euler_physics->compute_pressure(soln_at_iquad);
294  }
295 
296  // Obtain value used to linearly scale state variables
297  std::vector<real> theta2 = get_theta2_Zhang2010(p_lim, soln_cell_avg, soln_at_q, n_quad_pts, lower_bound, euler_physics->gam);
298 
299  // Limit values at quadrature points
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];
304  }
305  }
306  }
307 
308  else if (limiter_type == limiter_enum::positivity_preservingWang2012 && nstate == dim + 2) {
309  real theta2 = get_theta2_Wang2012(soln_at_q, n_quad_pts, p_avg); // Value used to linearly scale state variables
310 
311  // Limit values at quadrature points
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];
316  }
317  }
318  }
319 
320  // Project soln at quadrature points to dofs.
321  for (int istate = 0; istate < nstate; istate++) {
322  soln_basis_projection_oper.matrix_vector_mult_1D(soln_at_q[istate], soln_coeff[istate], soln_basis_projection_oper.oneD_vol_operator);
323  }
324 
325  // Write limited solution back and verify that positivity of density is satisfied
326  write_limited_solution(solution, soln_coeff, n_shape_fns, current_dofs_indices);
327  }
328 }
329 
336 } // PHiLiP namespace
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.
Definition: operators.cpp:1082
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 > &current_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.
Definition: operators.cpp:1732
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.
Definition: ADTypes.hpp:10
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.
Definition: operators.h:355
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.
Definition: operators.cpp:308
const Parameters::AllParameters *const all_parameters
Pointer to parameters object.
Projection operator corresponding to basis functions onto M-norm (L2).
Definition: operators.h:694