3 #include "convection_diffusion.h" 8 template <
int nstate,
typename real>
9 std::array<real,nstate> stdvector_to_stdarray(
const std::vector<real> vector)
11 std::array<real,nstate> array;
12 for (
int i=0; i<nstate; i++) { array[i] = vector[i]; }
16 template <
int dim,
int nstate,
typename real>
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 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);
34 for (
int istate=0; istate<nstate; ++istate) {
36 std::array<real,nstate> characteristic_dot_n = convective_eigenvalues(boundary_values, normal_int);
37 const bool inflow = (characteristic_dot_n[istate] <= 0.);
39 if (inflow || hasDiffusion) {
43 soln_bc[istate] = boundary_values[istate];
44 soln_grad_bc[istate] = soln_grad_int[istate];
50 soln_bc[istate] = soln_int[istate];
57 soln_grad_bc[istate] = soln_grad_int[istate];
64 template <
int dim,
int nstate,
typename real>
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) {
72 for (
int d=0; d<dim; ++d) {
73 conv_flux[i][d] += velocity_field[d] * solution[i];
79 template <
int dim,
int nstate,
typename real>
82 const std::array<real,nstate> &soln1,
83 const std::array<real,nstate> &soln2)
const 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;
89 return convective_flux(arr_avg);
92 template <
int dim,
int nstate,
typename real>
95 const std::array<real,nstate> &conservative_soln)
const 97 return conservative_soln;
100 template <
int dim,
int nstate,
typename real>
103 const std::array<real,nstate> &entropy_var)
const 108 template <
int dim,
int nstate,
typename real>
112 dealii::Tensor<1,dim,real> advection_speed;
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];
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;
123 return advection_speed;
126 template <
int dim,
int nstate,
typename real>
130 if(hasDiffusion)
return diffusion_scaling_coeff;
131 const real zero = 0.0;
135 template <
int dim,
int nstate,
typename real>
138 const std::array<real,nstate> &,
139 const dealii::Tensor<1,dim,real> &normal)
const 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];
147 for (
int i=0; i<nstate; i++) {
153 template <
int dim,
int nstate,
typename real>
157 const dealii::Tensor<1,dim,real> advection_speed = this->advection_speed();
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);
166 template <
int dim,
int nstate,
typename real>
170 const real diff_coeff = this->diffusion_coefficient();
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);
181 template <
int dim,
int nstate,
typename real>
184 const std::array<real,nstate> &,
185 const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient)
const 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]);
200 template <
int dim,
int nstate,
typename 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 )
const 207 return dissipative_flux(solution, solution_gradient);
210 template <
int dim,
int nstate,
typename 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 )
const 218 return source_term(pos,solution,current_time);
221 template <
int dim,
int nstate,
typename real>
224 const dealii::Point<dim,real> &pos,
225 const std::array<real,nstate> &,
226 const real current_time)
const 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();
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]);
243 source[istate] += (- diff_coeff) * exp(-diff_coeff * current_time) * sine_term;
244 for(
int idim=0; idim<dim; idim++){
245 source[istate] += diff_coeff * pow(pi,2) * exp(-diff_coeff * current_time)
246 * this->diffusion_tensor[idim][idim] * sine_term;
248 for(
int idim=0; idim<dim; idim++){
249 for(
int jdim=0; jdim<dim; jdim++){
251 real cross_term = cos(pi*pos[idim]) * cos(pi*pos[jdim]);
253 int kdim = 3 - idim - jdim;
254 cross_term *= sin(pi*pos[kdim]);
256 source[istate] += - diff_coeff * pow(pi,2) * exp(-diff_coeff * current_time)
257 * this->diffusion_tensor[idim][jdim] * cross_term;
264 for (
int istate=0; istate<nstate; istate++) {
265 dealii::Tensor<1,dim,real> manufactured_gradient = this->manufactured_solution_function->gradient (pos, istate);
273 dealii::SymmetricTensor<2,dim,real> manufactured_hessian = this->manufactured_solution_function->hessian (pos, istate);
284 for (
int d=0; d<dim; ++d) {
285 grad += velocity_field[d] * manufactured_gradient[d];
287 source[istate] = grad;
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];
295 source[istate] += -diff_coeff*hess;
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 'c'.
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.
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: .