[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
convection_diffusion.cpp
1 #include "ADTypes.hpp"
2 
3 #include "convection_diffusion.h"
4 
5 namespace PHiLiP {
6 namespace Physics {
7 
8 template <int nstate, typename real>
9 std::array<real,nstate> stdvector_to_stdarray(const std::vector<real> vector)
10 {
11  std::array<real,nstate> array;
12  for (int i=0; i<nstate; i++) { array[i] = vector[i]; }
13  return array;
14 }
15 
16 template <int dim, int nstate, typename real>
19  const int /*boundary_type*/,
20  const dealii::Point<dim, real> &pos,
21  const dealii::Tensor<1,dim,real> &normal_int,
22  const std::array<real,nstate> &soln_int,
23  const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
24  std::array<real,nstate> &soln_bc,
25  std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
26 {
27  std::array<real,nstate> boundary_values;
28  std::array<dealii::Tensor<1,dim,real>,nstate> boundary_gradients;
29  for (int i=0; i<nstate; i++) {
30  boundary_values[i] = this->manufactured_solution_function->value (pos, i);
31  boundary_gradients[i] = this->manufactured_solution_function->gradient (pos, i);
32  }
33 
34  for (int istate=0; istate<nstate; ++istate) {
35 
36  std::array<real,nstate> characteristic_dot_n = convective_eigenvalues(boundary_values, normal_int);
37  const bool inflow = (characteristic_dot_n[istate] <= 0.);
38 
39  if (inflow || hasDiffusion) { // Dirichlet boundary condition
40  // soln_bc[istate] = boundary_values[istate];
41  // soln_grad_bc[istate] = soln_grad_int[istate];
42 
43  soln_bc[istate] = boundary_values[istate];
44  soln_grad_bc[istate] = soln_grad_int[istate];
45 
46  } else { // Neumann boundary condition
47  // //soln_bc[istate] = soln_int[istate];
48  // //soln_bc[istate] = boundary_values[istate];
49  // soln_bc[istate] = -soln_int[istate]+2*boundary_values[istate];
50  soln_bc[istate] = soln_int[istate];
51 
52  // **************************************************************************************************************
53  // Note I don't know how to properly impose the soln_grad_bc to obtain an adjoint consistent scheme
54  // Currently, Neumann boundary conditions are only imposed for the linear advection
55  // Therefore, soln_grad_bc does not affect the solution
56  // **************************************************************************************************************
57  soln_grad_bc[istate] = soln_grad_int[istate];
58  //soln_grad_bc[istate] = boundary_gradients[istate];
59  //soln_grad_bc[istate] = -soln_grad_int[istate]+2*boundary_gradients[istate];
60  }
61  }
62 }
63 
64 template <int dim, int nstate, typename real>
65 std::array<dealii::Tensor<1,dim,real>,nstate> ConvectionDiffusion<dim,nstate,real>
66 ::convective_flux (const std::array<real,nstate> &solution) const
67 {
68  std::array<dealii::Tensor<1,dim,real>,nstate> conv_flux;
69  const dealii::Tensor<1,dim,real> velocity_field = advection_speed();
70  for (int i=0; i<nstate; ++i) {
71  conv_flux[i] = 0.0;
72  for (int d=0; d<dim; ++d) {
73  conv_flux[i][d] += velocity_field[d] * solution[i];
74  }
75  }
76  return conv_flux;
77 }
78 
79 template <int dim, int nstate, typename real>
80 std::array<dealii::Tensor<1,dim,real>,nstate> ConvectionDiffusion<dim,nstate,real>
82  const std::array<real,nstate> &soln1,
83  const std::array<real,nstate> &soln2) const
84 {
85  std::array<real,nstate> arr_avg;
86  for (int i = 0 ; i < nstate; ++i) {
87  arr_avg[i] = (soln1[i] + soln2[i])/2.0;
88  }
89  return convective_flux(arr_avg);
90 }
91 
92 template <int dim, int nstate, typename real>
93 std::array<real,nstate> ConvectionDiffusion<dim, nstate, real>
95  const std::array<real,nstate> &conservative_soln) const
96 {
97  return conservative_soln;
98 }
99 
100 template <int dim, int nstate, typename real>
101 std::array<real,nstate> ConvectionDiffusion<dim, nstate, real>
103  const std::array<real,nstate> &entropy_var) const
104 {
105  return entropy_var;
106 }
107 
108 template <int dim, int nstate, typename real>
109 dealii::Tensor<1,dim,real> ConvectionDiffusion<dim,nstate,real>
111 {
112  dealii::Tensor<1,dim,real> advection_speed;
113  if (hasConvection) {
114  if(dim >= 1) advection_speed[0] = linear_advection_velocity[0];
115  if(dim >= 2) advection_speed[1] = linear_advection_velocity[1];
116  if(dim >= 3) advection_speed[2] = linear_advection_velocity[2];
117  } else {
118  const real zero = 0.0;
119  if(dim >= 1) advection_speed[0] = zero;
120  if(dim >= 2) advection_speed[1] = zero;
121  if(dim >= 3) advection_speed[2] = zero;
122  }
123  return advection_speed;
124 }
125 
126 template <int dim, int nstate, typename real>
129 {
130  if(hasDiffusion) return diffusion_scaling_coeff;
131  const real zero = 0.0;
132  return zero;
133 }
134 
135 template <int dim, int nstate, typename real>
136 std::array<real,nstate> ConvectionDiffusion<dim,nstate,real>
138  const std::array<real,nstate> &/*solution*/,
139  const dealii::Tensor<1,dim,real> &normal) const
140 {
141  std::array<real,nstate> eig;
142  const dealii::Tensor<1,dim,real> advection_speed = this->advection_speed();
143  real eig_value = 0.0;
144  for (int d=0; d<dim; ++d) {
145  eig_value += advection_speed[d]*normal[d];
146  }
147  for (int i=0; i<nstate; i++) {
148  eig[i] = eig_value;
149  }
150  return eig;
151 }
152 
153 template <int dim, int nstate, typename real>
155 ::max_convective_eigenvalue (const std::array<real,nstate> &/*soln*/) const
156 {
157  const dealii::Tensor<1,dim,real> advection_speed = this->advection_speed();
158  real max_eig = 0;
159  for (int i=0; i<dim; i++) {
160  real abs_adv = abs(advection_speed[i]);
161  max_eig = std::max(max_eig,abs_adv);
162  }
163  return max_eig;
164 }
165 
166 template <int dim, int nstate, typename real>
168 ::max_viscous_eigenvalue (const std::array<real,nstate> &/*soln*/) const
169 {
170  const real diff_coeff = this->diffusion_coefficient();
171  real max_eig = 0;
172  for (int i=0; i<dim; i++) {
173  for (int j=0; j<dim; j++) {
174  real abs_visc = abs(diff_coeff * this->diffusion_tensor[i][j]);
175  max_eig = std::max(max_eig,abs_visc);
176  }
177  }
178  return max_eig;
179 }
180 
181 template <int dim, int nstate, typename real>
182 std::array<dealii::Tensor<1,dim,real>,nstate> ConvectionDiffusion<dim,nstate,real>
184  const std::array<real,nstate> &/*solution*/,
185  const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient) const
186 {
187  std::array<dealii::Tensor<1,dim,real>,nstate> diss_flux;
188  const real diff_coeff = diffusion_coefficient();
189  for (int i=0; i<nstate; i++) {
190  for (int d1=0; d1<dim; d1++) {
191  diss_flux[i][d1] = 0.0;
192  for (int d2=0; d2<dim; d2++) {
193  diss_flux[i][d1] += -diff_coeff*(this->diffusion_tensor[d1][d2]*solution_gradient[i][d2]);
194  }
195  }
196  }
197  return diss_flux;
198 }
199 
200 template <int dim, int nstate, typename real>
201 std::array<dealii::Tensor<1,dim,real>,nstate> ConvectionDiffusion<dim,nstate,real>
203  const std::array<real,nstate> &solution,
204  const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
205  const dealii::types::global_dof_index /*cell_index*/) const
206 {
207  return dissipative_flux(solution, solution_gradient);
208 }
209 
210 template <int dim, int nstate, typename real>
211 std::array<real,nstate> ConvectionDiffusion<dim,nstate,real>
213  const dealii::Point<dim,real> &pos,
214  const std::array<real,nstate> &solution,
215  const real current_time,
216  const dealii::types::global_dof_index /*cell_index*/) const
217 {
218  return source_term(pos,solution,current_time);
219 }
220 
221 template <int dim, int nstate, typename real>
222 std::array<real,nstate> ConvectionDiffusion<dim,nstate,real>
224  const dealii::Point<dim,real> &pos,
225  const std::array<real,nstate> &/*solution*/,
226  const real current_time) const
227 {
228  std::array<real,nstate> source;
229  const dealii::Tensor<1,dim,real> velocity_field = this->advection_speed();
230  const real diff_coeff = diffusion_coefficient();
231 
232 
233  using TestType = Parameters::AllParameters::TestType;
234 
235  if(this->test_type == TestType::convection_diffusion_periodicity){
236  for(int istate =0; istate<nstate; istate++){
237  source[istate] = 0.0;
238  const double pi = atan(1)*4.0;
239  real sine_term = 1.0;
240  for(int idim=0; idim<dim; idim++){
241  sine_term *= sin(pi * pos[idim]);
242  }
243  source[istate] += (- diff_coeff) * exp(-diff_coeff * current_time) * sine_term;//the unsteady term
244  for(int idim=0; idim<dim; idim++){//laplacian term
245  source[istate] += diff_coeff * pow(pi,2) * exp(-diff_coeff * current_time)
246  * this->diffusion_tensor[idim][idim] * sine_term;
247  }
248  for(int idim=0; idim<dim; idim++){//cross terms
249  for(int jdim=0; jdim<dim; jdim++){
250  if(idim != jdim){
251  real cross_term = cos(pi*pos[idim]) * cos(pi*pos[jdim]);
252  if(dim == 3){
253  int kdim = 3 - idim - jdim;
254  cross_term *= sin(pi*pos[kdim]);
255  }
256  source[istate] += - diff_coeff * pow(pi,2) * exp(-diff_coeff * current_time)
257  * this->diffusion_tensor[idim][jdim] * cross_term;
258  }
259  }
260  }
261  }
262  }
263  else{
264  for (int istate=0; istate<nstate; istate++) {
265  dealii::Tensor<1,dim,real> manufactured_gradient = this->manufactured_solution_function->gradient (pos, istate);
266  // dealii::Tensor<1,dim,real> manufactured_gradient_fd = this->manufactured_solution_function.gradient_fd (pos, istate);
267  // std::cout<<"FD" <<std::endl;
268  // std::cout<<manufactured_gradient_fd <<std::endl;
269  // std::cout<<"AN" <<std::endl;
270  // std::cout<<manufactured_gradient <<std::endl;
271  // std::cout<<"DIFF" <<std::endl;
272  // std::cout<<manufactured_gradient - manufactured_gradient_fd <<std::endl;
273  dealii::SymmetricTensor<2,dim,real> manufactured_hessian = this->manufactured_solution_function->hessian (pos, istate);
274  // dealii::SymmetricTensor<2,dim,real> manufactured_hessian_fd = this->manufactured_solution_function.hessian_fd (pos, istate);
275  // std::cout<<"FD" <<std::endl;
276  // std::cout<<manufactured_hessian_fd <<std::endl;
277  // std::cout<<"AN" <<std::endl;
278  // std::cout<<manufactured_hessian <<std::endl;
279  // std::cout<<"DIFF" <<std::endl;
280  // std::cout<<manufactured_hessian - manufactured_hessian_fd <<std::endl;
281 
282  //source[istate] = velocity_field*manufactured_gradient;
283  real grad = 0.0;
284  for (int d=0; d<dim; ++d) {
285  grad += velocity_field[d] * manufactured_gradient[d];
286  }
287  source[istate] = grad;
288 
289  real hess = 0.0;
290  for (int dr=0; dr<dim; ++dr) {
291  for (int dc=0; dc<dim; ++dc) {
292  hess += (this->diffusion_tensor)[dr][dc] * manufactured_hessian[dr][dc];
293  }
294  }
295  source[istate] += -diff_coeff*hess;
296  }
297  }
298  return source;
299 }
300 
307 
314 
321 
328 
335 
336 } // Physics namespace
337 } // PHiLiP namespace
338 
339 
std::array< real, nstate > convective_eigenvalues(const std::array< real, nstate > &, const dealii::Tensor< 1, dim, real > &) const
Spectral radius of convective term Jacobian is &#39;c&#39;.
std::array< real, nstate > compute_entropy_variables(const std::array< real, nstate > &conservative_soln) const
Computes the entropy variables.
TestType
Possible integration tests to run.
Convection-diffusion with linear advective and diffusive term. Derived from PhysicsBase.
std::array< dealii::Tensor< 1, dim, real >, nstate > convective_numerical_split_flux(const std::array< real, nstate > &soln1, const std::array< real, nstate > &soln2) const override
Convective numerical split flux for split form.
Files for the baseline physics.
Definition: ADTypes.hpp:10
void boundary_face_values(const int, const dealii::Point< dim, real > &, const dealii::Tensor< 1, dim, real > &, const std::array< real, nstate > &, const std::array< dealii::Tensor< 1, dim, real >, nstate > &, std::array< real, nstate > &, std::array< dealii::Tensor< 1, dim, real >, nstate > &) const
If diffusion is present, assign Dirichlet boundary condition.
std::array< real, nstate > source_term(const dealii::Point< dim, real > &pos, const std::array< real, nstate > &solution, const real current_time, const dealii::types::global_dof_index cell_index) const
Source term is zero or depends on manufactured solution.
real max_viscous_eigenvalue(const std::array< real, nstate > &soln) const
Maximum viscous eigenvalue.
std::array< real, nstate > compute_conservative_variables_from_entropy_variables(const std::array< real, nstate > &entropy_var) const
Computes the conservative variables from the entropy variables.
real max_convective_eigenvalue(const std::array< real, nstate > &soln) const
Maximum convective eigenvalue.
dealii::Tensor< 1, dim, real > advection_speed() const
Linear advection speed: c.
std::array< dealii::Tensor< 1, dim, real >, nstate > dissipative_flux(const std::array< real, nstate > &solution, const std::array< dealii::Tensor< 1, dim, real >, nstate > &solution_gradient, const dealii::types::global_dof_index cell_index) const
Dissipative flux: u.
real diffusion_coefficient() const
Diffusion coefficient.
std::array< dealii::Tensor< 1, dim, real >, nstate > convective_flux(const std::array< real, nstate > &solution) const
Convective flux: .