1 #include <CoDiPack/include/codi.hpp> 3 #include <deal.II/base/function.h> 4 #include <deal.II/base/function.templates.h> 5 #include <deal.II/base/function_time.templates.h> 7 #include "manufactured_solution.h" 9 template class dealii::FunctionTime<Sacado::Fad::DFad<double>>;
10 template class dealii::Function<PHILIP_DIM,Sacado::Fad::DFad<double>>;
17 return std::isfinite(static_cast<double>(value));
23 return std::isfinite(static_cast<double>(value.val()));
27 bool isfinite(Sacado::Fad::DFad<Sacado::Fad::DFad<double>> value)
29 return std::isfinite(static_cast<double>(value.val().val()));
33 bool isfinite(Sacado::Rad::ADvar<Sacado::Fad::DFad<double>> value)
35 return std::isfinite(static_cast<double>(value.val().val()));
38 template <
int dim,
typename real>
40 ::value (
const dealii::Point<dim,real> &,
const unsigned int )
const 46 template <
int dim,
typename real>
48 ::value (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 50 real value = this->amplitudes[istate];
51 for (
int d=0; d<dim; d++) {
52 value *= sin( this->frequencies[istate][d] * point[d] );
55 value += this->base_values[istate];
59 template <
int dim,
typename real>
61 ::value (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 64 for (
int d=0; d<dim; d++) {
65 value += this->amplitudes[istate]*sin( this->frequencies[istate][d] * point[d] );
68 value += this->base_values[istate];
72 template <
int dim,
typename real>
74 ::value (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 76 real value = this->amplitudes[istate];
77 for (
int d=0; d<dim; d++) {
78 value *= cos( this->frequencies[istate][d] * point[d] );
81 value += this->base_values[istate];
85 template <
int dim,
typename real>
87 ::value (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 90 for (
int d=0; d<dim; d++) {
91 value += exp( point[d] );
94 value += this->base_values[istate];
98 template <
int dim,
typename real>
100 ::value (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 103 const double poly_max = 7;
104 for (
int d=0; d<dim; d++) {
105 value += pow(point[d] + 0.5, poly_max);
107 value += this->base_values[istate];
111 template <
int dim,
typename real>
113 ::value (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 116 for (
int d=0; d<dim; d++) {
117 const real x = point[d];
118 value += 1.0 + x - x*x - x*x*x + x*x*x*x - x*x*x*x*x + x*x*x*x*x*x + 0.001*sin(50*x);
120 value += this->base_values[istate];
124 template <
int dim,
typename real>
126 ::value(
const dealii::Point<dim,real> &point,
const unsigned int )
const 129 for(
unsigned int i = 0; i < dim; ++i){
132 for(
unsigned int j = 0; j < n_shocks[i]; ++j){
134 val_dim += atan(S_j[i][j]*(x-x_j[i][j]));
141 template <
int dim,
typename real>
143 ::value(
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 146 for(
unsigned int d = 0; d < dim; ++d){
148 val *= x + (exp(x/epsilon[istate][d])-1.0)/(1.0-exp(1.0/epsilon[istate][d]));
153 template <
int dim,
typename real>
155 ::value(
const dealii::Point<dim,real> &point,
const unsigned int )
const 159 const real x = point[0], y = point[1];
162 val = a*tanh(b*sin(c*y + d) + e*x + f);
167 template <
int dim,
typename real>
169 ::value(
const dealii::Point<dim,real> &point,
const unsigned int )
const 173 for(
unsigned int d = 0; d < dim; ++d){
175 val += alpha_diag[d]*x*x;
180 template <
int dim,
typename real>
182 ::value (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 185 for (
int d=0; d<dim; d++) {
186 value += exp( point[d] ) + sin(point[d] );
189 value += this->base_values[istate];
193 template <
int dim,
typename real>
198 if constexpr(dim == 2) {
199 const real x = point[0], y = point[1];
205 value = ncm[0][0] + ncm[0][1]*sin(ncm[0][4]*c*x) + ncm[0][2]*cos(ncm[0][5]*c*y) + ncm[0][3]*cos(ncm[0][6]*c*x)*cos(ncm[0][6]*c*y);
209 value = ncm[1][0] + ncm[1][1]*sin(ncm[1][4]*c*x) + ncm[1][2]*cos(ncm[1][5]*c*y) + ncm[1][3]*cos(ncm[1][6]*c*x)*cos(ncm[1][6]*c*y);
213 value = ncm[2][0] + ncm[2][1]*cos(ncm[2][4]*c*x) + ncm[2][2]*sin(ncm[2][5]*c*y) + ncm[2][3]*cos(ncm[2][6]*c*x)*cos(ncm[2][6]*c*y);
217 value = ncm[3][0] + ncm[3][1]*cos(ncm[3][4]*c*x) + ncm[3][2]*sin(ncm[3][5]*c*y) + ncm[3][3]*cos(ncm[3][6]*c*x)*cos(ncm[3][6]*c*y);
221 value = ncm[4][0] + ncm[4][1]*cos(ncm[4][4]*c*x) + ncm[4][2]*cos(ncm[4][5]*c*y) + ncm[4][3]*cos(ncm[4][6]*c*x)*cos(ncm[4][6]*c*y);
227 template <
int dim,
typename real>
229 ::value(
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 233 const real density = primitive_value(point,0);
234 const real x_velocity = primitive_value(point,1);
235 const real y_velocity = primitive_value(point,2);
236 const real pressure = primitive_value(point,3);
237 const real turbulent_working_variable = primitive_value(point,4);
241 if(istate==0) value = density;
243 if(istate==1) value = density*x_velocity;
245 if(istate==2) value = density*y_velocity;
247 if(istate==3) value = pressure/(1.4-1.0) + 0.5*density*(x_velocity*x_velocity + y_velocity*y_velocity);
249 if(istate==4) value = density*turbulent_working_variable;
254 template <
int dim,
typename real>
256 ::gradient (
const dealii::Point<dim,real> &,
const unsigned int )
const 258 dealii::Tensor<1,dim,real> gradient;
259 for(
unsigned int i = 0; i < dim; i++){
265 template <
int dim,
typename real>
267 ::gradient (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 269 dealii::Tensor<1,dim,real> gradient;
270 for (
int dim_deri=0; dim_deri<dim; dim_deri++) {
271 gradient[dim_deri] = this->amplitudes[istate] * this->frequencies[istate][dim_deri];
272 for (
int dim_trig=0; dim_trig<dim; dim_trig++) {
273 const real angle = this->frequencies[istate][dim_trig] * point[dim_trig];
274 if (dim_deri == dim_trig) gradient[dim_deri] *= cos( angle );
275 if (dim_deri != dim_trig) gradient[dim_deri] *= sin( angle );
277 assert(
isfinite(gradient[dim_deri]));
280 const real A = this->amplitudes[istate];
281 const dealii::Tensor<1,dim,real> f = this->frequencies[istate];
283 const real fx = f[0]*point[0];
284 gradient[0] = A*f[0]*cos(fx);
287 const real fx = f[0]*point[0];
288 const real fy = f[1]*point[1];
289 gradient[0] = A*f[0]*cos(fx)*sin(fy);
290 gradient[1] = A*f[1]*sin(fx)*cos(fy);
293 const real fx = f[0]*point[0];
294 const real fy = f[1]*point[1];
295 const real fz = f[2]*point[2];
296 gradient[0] = A*f[0]*cos(fx)*sin(fy)*sin(fz);
297 gradient[1] = A*f[1]*sin(fx)*cos(fy)*sin(fz);
298 gradient[2] = A*f[2]*sin(fx)*sin(fy)*cos(fz);
303 template <
int dim,
typename real>
305 ::gradient (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 307 dealii::Tensor<1,dim,real> gradient;
308 const real A = this->amplitudes[istate];
309 const dealii::Tensor<1,dim,real> f = this->frequencies[istate];
311 const real fx = f[0]*point[0];
312 gradient[0] = A*f[0]*cos(fx);
315 const real fx = f[0]*point[0];
316 const real fy = f[1]*point[1];
317 gradient[0] = A*f[0]*cos(fx);
318 gradient[1] = A*f[1]*cos(fy);
321 const real fx = f[0]*point[0];
322 const real fy = f[1]*point[1];
323 const real fz = f[2]*point[2];
324 gradient[0] = A*f[0]*cos(fx);
325 gradient[1] = A*f[1]*cos(fy);
326 gradient[2] = A*f[2]*cos(fz);
331 template <
int dim,
typename real>
333 ::gradient (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 335 dealii::Tensor<1,dim,real> gradient;
336 const real A = this->amplitudes[istate];
337 const dealii::Tensor<1,dim,real> f = this->frequencies[istate];
339 const real fx = f[0]*point[0];
340 gradient[0] = -A*f[0]*sin(fx);
343 const real fx = f[0]*point[0];
344 const real fy = f[1]*point[1];
345 gradient[0] = -A*f[0]*sin(fx)*cos(fy);
346 gradient[1] = -A*f[1]*cos(fx)*sin(fy);
349 const real fx = f[0]*point[0];
350 const real fy = f[1]*point[1];
351 const real fz = f[2]*point[2];
352 gradient[0] = -A*f[0]*sin(fx)*cos(fy)*cos(fz);
353 gradient[1] = -A*f[1]*cos(fx)*sin(fy)*cos(fz);
354 gradient[2] = -A*f[2]*cos(fx)*cos(fy)*sin(fz);
359 template <
int dim,
typename real>
361 ::gradient (
const dealii::Point<dim,real> &point,
const unsigned int )
const 363 dealii::Tensor<1,dim,real> gradient;
365 gradient[0] = exp(point[0]);
368 gradient[0] = exp(point[0]);
369 gradient[1] = exp(point[1]);
372 gradient[0] = exp(point[0]);
373 gradient[1] = exp(point[1]);
374 gradient[2] = exp(point[2]);
379 template <
int dim,
typename real>
381 ::gradient (
const dealii::Point<dim,real> &point,
const unsigned int )
const 383 dealii::Tensor<1,dim,real> gradient;
384 const double poly_max = 7;
386 gradient[0] = poly_max*pow(point[0] + 0.5, poly_max-1);
389 gradient[0] = poly_max*pow(point[0] + 0.5, poly_max-1);
390 gradient[1] = poly_max*pow(point[1] + 0.5, poly_max-1);
393 gradient[0] = poly_max*pow(point[0] + 0.5, poly_max-1);
394 gradient[1] = poly_max*pow(point[1] + 0.5, poly_max-1);
395 gradient[2] = poly_max*pow(point[2] + 0.5, poly_max-1);
400 template <
int dim,
typename real>
402 ::gradient (
const dealii::Point<dim,real> &point,
const unsigned int )
const 404 dealii::Tensor<1,dim,real> gradient;
406 const real x = point[0];
407 gradient[0] = 1.0 - 2*x -3*x*x + 4*x*x*x - 5*x*x*x*x + 6*x*x*x*x*x + 0.050*cos(50*x);
411 gradient[0] = 1.0 - 2*x -3*x*x + 4*x*x*x - 5*x*x*x*x + 6*x*x*x*x*x + 0.050*cos(50*x);
413 gradient[1] = 1.0 - 2*x -3*x*x + 4*x*x*x - 5*x*x*x*x + 6*x*x*x*x*x + 0.050*cos(50*x);
417 gradient[0] = 1.0 - 2*x -3*x*x + 4*x*x*x - 5*x*x*x*x + 6*x*x*x*x*x;
419 gradient[1] = 1.0 - 2*x -3*x*x + 4*x*x*x - 5*x*x*x*x + 6*x*x*x*x*x;
421 gradient[2] = 1.0 - 2*x -3*x*x + 4*x*x*x - 5*x*x*x*x + 6*x*x*x*x*x;
426 template <
int dim,
typename real>
428 ::gradient(
const dealii::Point<dim,real> &point,
const unsigned int )
const 430 dealii::Tensor<1,dim,real> gradient;
431 for(
unsigned int k = 0; k < dim; ++k){
434 for(
unsigned int i = 0; i < dim; ++i){
437 for(
unsigned int j = 0; j < n_shocks[i]; ++j){
440 real coeff = S_j[i][j]*(x-x_j[i][j]);
441 val_dim += S_j[i][j]/(pow(coeff,2)+1);
444 val_dim += atan(S_j[i][j]*(x-x_j[i][j]));
449 gradient[k] = grad_dim;
454 template <
int dim,
typename real>
456 ::gradient(
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 458 dealii::Tensor<1,dim,real> gradient;
460 const real x = point[0];
461 gradient[0] = (1 + (exp(x/epsilon[istate][0])/epsilon[istate][0])/(1.0-exp(1.0/epsilon[istate][0])));
463 const real x = point[0], y = point[1];
464 gradient[0] = (1 + (exp(x/epsilon[istate][0])/epsilon[istate][0])/(1.0-exp(1.0/epsilon[istate][0])))
465 * (y + (exp(y/epsilon[istate][1])-1.0) /(1.0-exp(1.0/epsilon[istate][1])));
466 gradient[1] = (x + (exp(x/epsilon[istate][0])-1.0) /(1.0-exp(1.0/epsilon[istate][0])))
467 * (1 + (exp(y/epsilon[istate][1])/epsilon[istate][1])/(1.0-exp(1.0/epsilon[istate][1])));
469 const real x = point[0], y = point[1], z = point[2];
470 gradient[0] = (1 + (exp(x/epsilon[istate][0])/epsilon[istate][0])/(1.0-exp(1.0/epsilon[istate][0])))
471 * (y + (exp(y/epsilon[istate][1])-1.0) /(1.0-exp(1.0/epsilon[istate][1])))
472 * (z + (exp(z/epsilon[istate][2])-1.0) /(1.0-exp(1.0/epsilon[istate][2])));
473 gradient[1] = (x + (exp(x/epsilon[istate][0])-1.0) /(1.0-exp(1.0/epsilon[istate][0])))
474 * (1 + (exp(y/epsilon[istate][1])/epsilon[istate][1])/(1.0-exp(1.0/epsilon[istate][1])))
475 * (z + (exp(z/epsilon[istate][2])-1.0) /(1.0-exp(1.0/epsilon[istate][2])));
476 gradient[2] = (x + (exp(x/epsilon[istate][0])-1.0) /(1.0-exp(1.0/epsilon[istate][0])))
477 * (y + (exp(y/epsilon[istate][1])-1.0) /(1.0-exp(1.0/epsilon[istate][1])))
478 * (1 + (exp(z/epsilon[istate][2])/epsilon[istate][2])/(1.0-exp(1.0/epsilon[istate][2])));
483 template <
int dim,
typename real>
485 ::gradient(
const dealii::Point<dim,real> &point,
const unsigned int )
const 487 dealii::Tensor<1,dim,real> gradient;
489 const real x = point[0], y = point[1];
495 const real denominator = pow(cosh(f + e*x + b*sin(d + c*y)), -2);
496 gradient[0] = a*e*denominator;
497 gradient[1] = a*b*c*cos(d+c*y)*denominator;
502 template <
int dim,
typename real>
504 ::gradient(
const dealii::Point<dim,real> &point,
const unsigned int )
const 506 dealii::Tensor<1,dim,real> gradient;
507 for(
unsigned int d = 0; d < dim; ++d){
509 const real x = point[d];
510 gradient[d] = 2*alpha_diag[d]*x;
515 template <
int dim,
typename real>
517 ::gradient(
const dealii::Point<dim,real> &point,
const unsigned int )
const 519 dealii::Tensor<1,dim,real> gradient;
520 for(
unsigned int d = 0; d < dim; ++d){
521 gradient[d] = exp(point[d]) + cos(point[d]);
526 template <
int dim,
typename real>
530 dealii::Tensor<1,dim,real> gradient;
533 const real x = point[0], y = point[1];
537 gradient[0] = ncm[0][4]*c*ncm[0][1]*cos(ncm[0][4]*c*x) - ncm[0][6]*c*ncm[0][3]*sin(ncm[0][6]*c*x)*cos(ncm[0][6]*c*y);
538 gradient[1] = -ncm[0][5]*c*ncm[0][2]*sin(ncm[0][5]*c*y) - ncm[0][6]*c*ncm[0][3]*cos(ncm[0][6]*c*x)*sin(ncm[0][6]*c*y);
542 gradient[0] = ncm[1][4]*c*ncm[1][1]*cos(ncm[1][4]*c*x) - ncm[1][6]*c*ncm[1][3]*sin(ncm[1][6]*c*x)*cos(ncm[1][6]*c*y);
543 gradient[1] = -ncm[1][5]*c*ncm[1][2]*sin(ncm[1][5]*c*y) - ncm[1][6]*c*ncm[1][3]*cos(ncm[1][6]*c*x)*sin(ncm[1][6]*c*y);
547 gradient[0] = -ncm[2][4]*c*ncm[2][1]*sin(ncm[2][4]*c*x) - ncm[2][6]*c*ncm[2][3]*sin(ncm[2][6]*c*x)*cos(ncm[2][6]*c*y);
548 gradient[1] = ncm[2][5]*c*ncm[2][2]*cos(ncm[2][5]*c*y) - ncm[2][6]*c*ncm[2][3]*cos(ncm[2][6]*c*x)*sin(ncm[2][6]*c*y);
552 gradient[0] = -ncm[3][4]*c*ncm[3][1]*sin(ncm[3][4]*c*x) - ncm[3][6]*c*ncm[3][3]*sin(ncm[3][6]*c*x)*cos(ncm[3][6]*c*y);
553 gradient[1] = ncm[3][5]*c*ncm[3][2]*cos(ncm[3][5]*c*y) - ncm[3][6]*c*ncm[3][3]*cos(ncm[3][6]*c*x)*sin(ncm[3][6]*c*y);
557 gradient[0] = -ncm[4][4]*c*ncm[4][1]*sin(ncm[4][4]*c*x) - ncm[4][6]*c*ncm[4][3]*sin(ncm[4][6]*c*x)*cos(ncm[4][6]*c*y);
558 gradient[1] = -ncm[4][5]*c*ncm[4][2]*sin(ncm[4][5]*c*y) - ncm[4][6]*c*ncm[4][3]*cos(ncm[4][6]*c*x)*sin(ncm[4][6]*c*y);
564 template <
int dim,
typename real>
566 ::gradient (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 568 dealii::Tensor<1,dim,real> gradient;
571 const real rho = primitive_value(point,0);
572 const real u = primitive_value(point,1);
573 const real v = primitive_value(point,2);
575 const dealii::Tensor<1,dim,real> rho_grad = primitive_gradient(point,0);
576 const dealii::Tensor<1,dim,real> u_grad = primitive_gradient(point,1);
577 const dealii::Tensor<1,dim,real> v_grad = primitive_gradient(point,2);
578 const dealii::Tensor<1,dim,real> p_grad = primitive_gradient(point,3);
583 for(
int d=0; d<dim; d++) {
584 gradient[d] = rho_grad[d];
589 for(
int d=0; d<dim; d++) {
590 gradient[d] = u*rho_grad[d] + rho*u_grad[d];
595 for(
int d=0; d<dim; d++) {
596 gradient[d] = v*rho_grad[d] + rho*v_grad[d];
601 for(
int d=0; d<dim; d++) {
602 gradient[d] = p_grad[d]/(1.4-1.0) + 0.5*rho_grad[d]*(u*u + v*v) + rho*(u*u_grad[d]+v*v_grad[d]);
607 const real twv = primitive_value(point,4);
608 const dealii::Tensor<1,dim,real> twv_grad = primitive_gradient(point,4);
609 for(
int d=0; d<dim; d++) {
610 gradient[d] = twv*rho_grad[d] + rho*twv_grad[d];
617 template <
int dim,
typename real>
619 ::hessian (
const dealii::Point<dim,real> &,
const unsigned int )
const 621 dealii::SymmetricTensor<2,dim,real> hessian;
622 for(
unsigned int i = 0; i < dim; i++){
623 for(
unsigned int j = 0; j < dim; j++){
630 template <
int dim,
typename real>
632 ::hessian (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 634 dealii::SymmetricTensor<2,dim,real> hessian;
636 const real A = this->amplitudes[istate];
637 const dealii::Tensor<1,dim,real> f = this->frequencies[istate];
639 const real fx = f[0]*point[0];
640 hessian[0][0] = -A*f[0]*f[0]*sin(fx);
643 const real fx = f[0]*point[0];
644 const real fy = f[1]*point[1];
645 hessian[0][0] = -A*f[0]*f[0]*sin(fx)*sin(fy);
646 hessian[0][1] = A*f[0]*f[1]*cos(fx)*cos(fy);
648 hessian[1][0] = A*f[1]*f[0]*cos(fx)*cos(fy);
649 hessian[1][1] = -A*f[1]*f[1]*sin(fx)*sin(fy);
652 const real fx = f[0]*point[0];
653 const real fy = f[1]*point[1];
654 const real fz = f[2]*point[2];
655 hessian[0][0] = -A*f[0]*f[0]*sin(fx)*sin(fy)*sin(fz);
656 hessian[0][1] = A*f[0]*f[1]*cos(fx)*cos(fy)*sin(fz);
657 hessian[0][2] = A*f[0]*f[2]*cos(fx)*sin(fy)*cos(fz);
659 hessian[1][0] = A*f[1]*f[0]*cos(fx)*cos(fy)*sin(fz);
660 hessian[1][1] = -A*f[1]*f[1]*sin(fx)*sin(fy)*sin(fz);
661 hessian[1][2] = A*f[1]*f[2]*sin(fx)*cos(fy)*cos(fz);
663 hessian[2][0] = A*f[2]*f[0]*cos(fx)*sin(fy)*cos(fz);
664 hessian[2][1] = A*f[2]*f[1]*sin(fx)*cos(fy)*cos(fz);
665 hessian[2][2] = -A*f[2]*f[2]*sin(fx)*sin(fy)*sin(fz);
670 template <
int dim,
typename real>
672 ::hessian (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 674 dealii::SymmetricTensor<2,dim,real> hessian;
675 const real A = this->amplitudes[istate];
676 const dealii::Tensor<1,dim,real> f = this->frequencies[istate];
678 const real fx = f[0]*point[0];
679 hessian[0][0] = -A*f[0]*f[0]*sin(fx);
682 const real fx = f[0]*point[0];
683 const real fy = f[1]*point[1];
684 hessian[0][0] = -A*f[0]*f[0]*sin(fx);
688 hessian[1][1] = -A*f[1]*f[1]*sin(fy);
691 const real fx = f[0]*point[0];
692 const real fy = f[1]*point[1];
693 const real fz = f[2]*point[2];
694 hessian[0][0] = -A*f[0]*f[0]*sin(fx);
699 hessian[1][1] = -A*f[1]*f[1]*sin(fy);
704 hessian[2][2] = -A*f[2]*f[2]*sin(fz);
709 template <
int dim,
typename real>
711 ::hessian (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 713 dealii::SymmetricTensor<2,dim,real> hessian;
714 const real A = this->amplitudes[istate];
715 const dealii::Tensor<1,dim,real> f = this->frequencies[istate];
717 const real fx = f[0]*point[0];
718 hessian[0][0] = -A*f[0]*f[0]*cos(fx);
721 const real fx = f[0]*point[0];
722 const real fy = f[1]*point[1];
723 hessian[0][0] = -A*f[0]*f[0]*cos(fx)*cos(fy);
724 hessian[0][1] = A*f[0]*f[1]*sin(fx)*sin(fy);
726 hessian[1][0] = A*f[1]*f[0]*sin(fx)*sin(fy);
727 hessian[1][1] = -A*f[1]*f[1]*cos(fx)*cos(fy);
730 const real fx = f[0]*point[0];
731 const real fy = f[1]*point[1];
732 const real fz = f[2]*point[2];
733 hessian[0][0] = -A*f[0]*f[0]*cos(fx)*cos(fy)*cos(fz);
734 hessian[0][1] = A*f[0]*f[1]*sin(fx)*sin(fy)*cos(fz);
735 hessian[0][2] = A*f[0]*f[2]*sin(fx)*cos(fy)*sin(fz);
737 hessian[1][0] = A*f[1]*f[0]*sin(fx)*sin(fy)*cos(fz);
738 hessian[1][1] = -A*f[1]*f[1]*cos(fx)*cos(fy)*cos(fz);
739 hessian[1][2] = A*f[1]*f[2]*cos(fx)*sin(fy)*sin(fz);
741 hessian[2][0] = A*f[2]*f[0]*sin(fx)*cos(fy)*sin(fz);
742 hessian[2][1] = A*f[2]*f[1]*cos(fx)*sin(fy)*sin(fz);
743 hessian[2][2] = -A*f[2]*f[2]*cos(fx)*cos(fy)*cos(fz);
748 template <
int dim,
typename real>
750 ::hessian (
const dealii::Point<dim,real> &point,
const unsigned int )
const 752 dealii::SymmetricTensor<2,dim,real> hessian;
754 hessian[0][0] = exp(point[0]);
757 hessian[0][0] = exp(point[0]);
761 hessian[1][1] = exp(point[1]);
764 hessian[0][0] = exp(point[0]);
769 hessian[1][1] = exp(point[1]);
774 hessian[2][2] = exp(point[2]);
779 template <
int dim,
typename real>
781 ::hessian (
const dealii::Point<dim,real> &point,
const unsigned int )
const 783 dealii::SymmetricTensor<2,dim,real> hessian;
784 const double poly_max = 7;
786 hessian[0][0] = poly_max*poly_max*pow(point[0] + 0.5, poly_max-2);
789 hessian[0][0] = poly_max*poly_max*pow(point[0] + 0.5, poly_max-2);
793 hessian[1][1] = poly_max*poly_max*pow(point[1] + 0.5, poly_max-2);
796 hessian[0][0] = poly_max*poly_max*pow(point[0] + 0.5, poly_max-2);
801 hessian[1][1] = poly_max*poly_max*pow(point[1] + 0.5, poly_max-2);
806 hessian[2][2] = poly_max*poly_max*pow(point[2] + 0.5, poly_max-2);
811 template <
int dim,
typename real>
813 ::hessian (
const dealii::Point<dim,real> &point,
const unsigned int )
const 815 dealii::SymmetricTensor<2,dim,real> hessian;
817 const real x = point[0];
818 hessian[0][0] = - 2.0 -6*x + 12*x*x - 20*x*x*x + 30*x*x*x*x - 2.500*sin(50*x);
822 hessian[0][0] = - 2.0 -6*x + 12*x*x - 20*x*x*x + 30*x*x*x*x - 2.500*sin(50*x);
824 hessian[1][1] = - 2.0 -6*x + 12*x*x - 20*x*x*x + 30*x*x*x*x - 2.500*sin(50*x);
828 hessian[0][0] = - 2.0 -6*x + 12*x*x - 20*x*x*x + 30*x*x*x*x - 2.500*sin(50*x);
830 hessian[1][1] = - 2.0 -6*x + 12*x*x - 20*x*x*x + 30*x*x*x*x - 2.500*sin(50*x);
832 hessian[2][2] = - 2.0 -6*x + 12*x*x - 20*x*x*x + 30*x*x*x*x - 2.500*sin(50*x);
837 template <
int dim,
typename real>
839 ::hessian(
const dealii::Point<dim,real> &point,
const unsigned int )
const 841 dealii::SymmetricTensor<2,dim,real> hes;
843 for(
unsigned int k1 = 0; k1 < dim; ++k1){
845 for(
unsigned int k2 = 0; k2 < dim; ++k2){
848 for(
unsigned int i = 0; i < dim; ++i){
851 for(
unsigned int j = 0; j < n_shocks[i]; ++j){
852 if(i == k1 && i == k2){
854 real coeff = S_j[i][j]*(x-x_j[i][j]);
855 val_dim += -2.0*pow(S_j[i][j],2)*coeff/pow(pow(coeff,2)+1,2);
856 }
else if(i == k1 || i == k2){
858 real coeff = S_j[i][j]*(x-x_j[i][j]);
859 val_dim += S_j[i][j]/(pow(coeff,2)+1);
862 val_dim += atan(S_j[i][j]*(x-x_j[i][j]));
867 hes[k1][k2] = hes_dim;
874 template <
int dim,
typename real>
876 ::hessian(
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 878 dealii::SymmetricTensor<2,dim,real> hessian;
880 const real x = point[0];
881 hessian[0][0] = (exp(x/epsilon[istate][0])/pow(epsilon[istate][0],2)/(1.0-exp(1.0/epsilon[istate][0])));
884 real x = point[0], y = point[1];
885 hessian[0][0] = (exp(x/epsilon[istate][0])/pow(epsilon[istate][0],2)/(1.0-exp(1.0/epsilon[istate][0])))
886 * (y + (exp(y/epsilon[istate][1])-1.0) /(1.0-exp(1.0/epsilon[istate][1])));
887 hessian[0][1] = (1 + (exp(x/epsilon[istate][0])/epsilon[istate][0]) /(1.0-exp(1.0/epsilon[istate][0])))
888 * (1 + (exp(y/epsilon[istate][1])/epsilon[istate][1]) /(1.0-exp(1.0/epsilon[istate][1])));
890 hessian[1][0] = hessian[0][1];
891 hessian[1][1] = (x + (exp(x/epsilon[istate][0])-1.0) /(1.0-exp(1.0/epsilon[istate][0])))
892 * (exp(y/epsilon[istate][1])/pow(epsilon[istate][1],2)/(1.0-exp(1.0/epsilon[istate][1])));
895 real x = point[0], y = point[1], z = point[2];
896 hessian[0][0] = (exp(x/epsilon[istate][0])/pow(epsilon[istate][0],2)/(1.0-exp(1.0/epsilon[istate][0])))
897 * (y + (exp(y/epsilon[istate][1])-1.0) /(1.0-exp(1.0/epsilon[istate][1])))
898 * (z + (exp(z/epsilon[istate][2])-1.0) /(1.0-exp(1.0/epsilon[istate][2])));
899 hessian[0][1] = (1 + (exp(x/epsilon[istate][0])/epsilon[istate][0]) /(1.0-exp(1.0/epsilon[istate][0])))
900 * (1 + (exp(y/epsilon[istate][1])/epsilon[istate][1]) /(1.0-exp(1.0/epsilon[istate][1])))
901 * (z + (exp(z/epsilon[istate][2])-1.0) /(1.0-exp(1.0/epsilon[istate][2])));
902 hessian[0][2] = (1 + (exp(x/epsilon[istate][0])/epsilon[istate][0]) /(1.0-exp(1.0/epsilon[istate][0])))
903 * (y + (exp(y/epsilon[istate][1])-1.0) /(1.0-exp(1.0/epsilon[istate][1])))
904 * (1 + (exp(z/epsilon[istate][2])/epsilon[istate][2]) /(1.0-exp(1.0/epsilon[istate][2])));
906 hessian[1][0] = hessian[0][1];
907 hessian[1][1] = (x + (exp(x/epsilon[istate][0])-1.0) /(1.0-exp(1.0/epsilon[istate][0])))
908 * (exp(y/epsilon[istate][1])/pow(epsilon[istate][1],2)/(1.0-exp(1.0/epsilon[istate][1])))
909 * (z + (exp(z/epsilon[istate][2])-1.0) /(1.0-exp(1.0/epsilon[istate][2])));
910 hessian[1][2] = (x + (exp(x/epsilon[istate][0])-1.0) /(1.0-exp(1.0/epsilon[istate][0])))
911 * (1 + (exp(y/epsilon[istate][1])/epsilon[istate][1]) /(1.0-exp(1.0/epsilon[istate][1])))
912 * (1 + (exp(z/epsilon[istate][2])/epsilon[istate][2]) /(1.0-exp(1.0/epsilon[istate][2])));
914 hessian[2][0] = hessian[0][2];
915 hessian[2][1] = hessian[2][1];
916 hessian[2][2] = (x + (exp(x/epsilon[istate][0])-1.0) /(1.0-exp(1.0/epsilon[istate][0])))
917 * (y + (exp(y/epsilon[istate][1])-1.0) /(1.0-exp(1.0/epsilon[istate][1])))
918 * (exp(z/epsilon[istate][2])/pow(epsilon[istate][2],2)/(1.0-exp(1.0/epsilon[istate][2])));
923 template <
int dim,
typename real>
925 ::hessian(
const dealii::Point<dim,real> &point,
const unsigned int )
const 927 dealii::SymmetricTensor<2,dim,real> hessian;
929 const real x = point[0], y = point[1];
942 const real component = f + e*x + b*sin(d+c*y);
943 const real numerator = sinh(component);
944 const real denominator = pow(cosh(component), -3);
946 hessian[0][0] = -2*a*e*e*numerator*denominator;
947 hessian[0][1] = -2*a*b*c*e*cos(d+c*y)*numerator*denominator;
949 hessian[1][0] = hessian[0][1];
950 hessian[1][1] = -a*b*c*c*pow(cosh(component), -2)*(2*b*pow(cos(c*y + d),2)*tanh(component) + sin(c*y + d));
955 template <
int dim,
typename real>
957 ::hessian(
const dealii::Point<dim,real> &,
const unsigned int )
const 959 dealii::SymmetricTensor<2,dim,real> hessian;
960 for(
unsigned int i = 0; i < dim; ++i){
961 for(
unsigned int j = 0; j < dim; ++j){
963 hessian[i][i] = 2*alpha_diag[i];
972 template <
int dim,
typename real>
974 ::hessian (
const dealii::Point<dim,real> &point,
const unsigned int )
const 976 dealii::SymmetricTensor<2,dim,real> hessian;
977 for(
int idim=0; idim<dim; idim++){
978 for(
int jdim=0; jdim<dim; jdim++){
980 hessian[idim][jdim] = exp(point[idim]) - sin(point[idim]);
982 hessian[idim][jdim] = 0.0;
988 template <
int dim,
typename real>
992 dealii::SymmetricTensor<2,dim,real> hessian;
995 const real x = point[0], y = point[1];
999 hessian[0][0] = -ncm[0][4]*c*ncm[0][4]*c*ncm[0][1]*sin(ncm[0][4]*c*x) - ncm[0][6]*c*ncm[0][6]*c*ncm[0][3]*cos(ncm[0][6]*c*x)*cos(ncm[0][6]*c*y);
1000 hessian[0][1] = ncm[0][6]*c*ncm[0][6]*c*ncm[0][3]*sin(ncm[0][6]*c*x)*sin(ncm[0][6]*c*y);
1001 hessian[1][0] = hessian[0][1];
1002 hessian[1][1] = -ncm[0][5]*c*ncm[0][5]*c*ncm[0][2]*cos(ncm[0][5]*c*y) - ncm[0][6]*c*ncm[0][6]*c*ncm[0][3]*cos(ncm[0][6]*c*x)*cos(ncm[0][6]*c*y);
1006 hessian[0][0] = -ncm[1][4]*c*ncm[1][4]*c*ncm[1][1]*sin(ncm[1][4]*c*x) - ncm[1][6]*c*ncm[1][6]*c*ncm[1][3]*cos(ncm[1][6]*c*x)*cos(ncm[1][6]*c*y);
1007 hessian[0][1] = ncm[1][6]*c*ncm[1][6]*c*ncm[1][3]*sin(ncm[1][6]*c*x)*sin(ncm[1][6]*c*y);
1008 hessian[1][0] = hessian[0][1];
1009 hessian[1][1] = -ncm[1][5]*c*ncm[1][5]*c*ncm[1][2]*cos(ncm[1][5]*c*y) - ncm[1][6]*c*ncm[1][6]*c*ncm[1][3]*cos(ncm[1][6]*c*x)*cos(ncm[1][6]*c*y);
1013 hessian[0][0] = -ncm[2][4]*c*ncm[2][4]*c*ncm[2][1]*cos(ncm[2][4]*c*x) - ncm[2][6]*c*ncm[2][6]*c*ncm[2][3]*cos(ncm[2][6]*c*x)*cos(ncm[2][6]*c*y);
1014 hessian[0][1] = ncm[2][6]*c*ncm[2][6]*c*ncm[2][3]*sin(ncm[2][6]*c*x)*sin(ncm[2][6]*c*y);
1015 hessian[1][0] = hessian[0][1];
1016 hessian[1][1] = -ncm[2][5]*c*ncm[2][5]*c*ncm[2][2]*sin(ncm[2][5]*c*y) - ncm[2][6]*c*ncm[2][6]*c*ncm[2][3]*cos(ncm[2][6]*c*x)*cos(ncm[2][6]*c*y);
1020 hessian[0][0] = -ncm[3][4]*c*ncm[3][4]*c*ncm[3][1]*cos(ncm[3][4]*c*x) - ncm[3][6]*c*ncm[3][6]*c*ncm[3][3]*cos(ncm[3][6]*c*x)*cos(ncm[3][6]*c*y);
1021 hessian[0][1] = ncm[3][6]*c*ncm[3][6]*c*ncm[3][3]*sin(ncm[3][6]*c*x)*sin(ncm[3][6]*c*y);
1022 hessian[1][0] = hessian[0][1];
1023 hessian[1][1] = -ncm[3][5]*c*ncm[3][5]*c*ncm[3][2]*sin(ncm[3][5]*c*y) - ncm[3][6]*c*ncm[3][6]*c*ncm[3][3]*cos(ncm[3][6]*c*x)*cos(ncm[3][6]*c*y);
1027 hessian[0][0] = -ncm[4][4]*c*ncm[4][4]*c*ncm[4][1]*cos(ncm[4][4]*c*x) - ncm[4][6]*c*ncm[4][6]*c*ncm[4][3]*cos(ncm[4][6]*c*x)*cos(ncm[4][6]*c*y);
1028 hessian[0][1] = ncm[4][6]*c*ncm[4][6]*c*ncm[4][3]*sin(ncm[4][6]*c*x)*sin(ncm[4][6]*c*y);
1029 hessian[1][0] = hessian[0][1];
1030 hessian[1][1] = -ncm[4][5]*c*ncm[4][5]*c*ncm[4][2]*cos(ncm[4][5]*c*y) - ncm[4][6]*c*ncm[4][6]*c*ncm[4][3]*cos(ncm[4][6]*c*x)*cos(ncm[4][6]*c*y);
1036 template <
int dim,
typename real>
1038 ::hessian (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 1040 dealii::SymmetricTensor<2,dim,real> hessian;
1043 const real rho = primitive_value(point,0);
1044 const real u = primitive_value(point,1);
1045 const real v = primitive_value(point,2);
1047 const dealii::Tensor<1,dim,real> rho_grad = primitive_gradient(point,0);
1048 const dealii::Tensor<1,dim,real> u_grad = primitive_gradient(point,1);
1049 const dealii::Tensor<1,dim,real> v_grad = primitive_gradient(point,2);
1051 const dealii::SymmetricTensor<2,dim,real> rho_hess = primitive_hessian(point,0);
1052 const dealii::SymmetricTensor<2,dim,real> u_hess = primitive_hessian(point,1);
1053 const dealii::SymmetricTensor<2,dim,real> v_hess = primitive_hessian(point,2);
1054 const dealii::SymmetricTensor<2,dim,real> p_hess = primitive_hessian(point,3);
1059 for(
int i=0; i<dim; i++) {
1060 for(
int j=0; j<dim; j++) {
1061 hessian[i][j] = rho_hess[i][j];
1067 for(
int i=0; i<dim; i++) {
1068 for(
int j=0; j<dim; j++) {
1069 hessian[i][j] = u_grad[j]*rho_grad[i] + u*rho_hess[i][j] + rho_grad[j]*u_grad[i] + rho*u_hess[i][j];
1075 for(
int i=0; i<dim; i++) {
1076 for(
int j=0; j<dim; j++) {
1077 hessian[i][j] = v_grad[j]*rho_grad[i] + v*rho_hess[i][j] + rho_grad[j]*v_grad[i] + rho*v_hess[i][j];
1083 for(
int i=0; i<dim; i++) {
1084 for(
int j=0; j<dim; j++) {
1085 hessian[i][j] = p_hess[i][j]/(1.4-1.0) + (u*u_grad[j]+v*v_grad[j])*rho_grad[i] + 0.5*(u*u + v*v)*rho_hess[i][j];
1086 hessian[i][j] += rho_grad[j]*(u*u_grad[i]+v*v_grad[i]);
1087 hessian[i][j] += rho*(u_grad[j]*u_grad[i] + u*u_hess[i][j] + v_grad[j]*v_grad[i] + v*v_hess[i][j]);
1093 const real twv = primitive_value(point,4);
1094 const dealii::Tensor<1,dim,real> twv_grad = primitive_gradient(point,4);
1095 const dealii::SymmetricTensor<2,dim,real> twv_hess = primitive_hessian(point,4);
1096 for(
int i=0; i<dim; i++) {
1097 for(
int j=0; j<dim; j++) {
1098 hessian[i][j] = twv_grad[j]*rho_grad[i] + twv*rho_hess[i][j] + rho_grad[j]*twv_grad[i] + rho*twv_hess[i][j];
1106 template <
int dim,
typename real>
1110 dealii::Function<dim,real>(nstate)
1112 , base_values(nstate)
1113 , amplitudes(nstate)
1114 , frequencies(nstate)
1116 const double pi = atan(1)*4.0;
1119 for (
int s=0; s<(int)nstate; s++) {
1120 base_values[s] = 1+(s+1.0)/nstate;
1121 base_values[nstate-1] = 10;
1122 amplitudes[s] = 0.2*base_values[s]*sin((static_cast<double>(nstate)-s)/nstate);
1123 for (
int d=0; d<dim; d++) {
1125 frequencies[s][d] = 2.0 + sin(0.1+s*0.5+d*0.2) * pi / 2.0;
1131 template <
int dim,
typename real>
1133 ::gradient_fd (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 1135 dealii::Tensor<1,dim,real> gradient;
1136 const double eps=1e-6;
1137 for (
int dim_deri=0; dim_deri<dim; dim_deri++) {
1138 dealii::Point<dim,real> pert_p = point;
1139 dealii::Point<dim,real> pert_m = point;
1140 pert_p[dim_deri] += eps;
1141 pert_m[dim_deri] -= eps;
1142 const real value_p = value(pert_p,istate);
1143 const real value_m = value(pert_m,istate);
1144 gradient[dim_deri] = (value_p - value_m) / (2*eps);
1149 template <
int dim,
typename real>
1151 ::hessian_fd (
const dealii::Point<dim,real> &point,
const unsigned int istate)
const 1153 dealii::SymmetricTensor<2,dim,real> hessian;
1154 const double eps=1e-4;
1155 for (
int d1=0; d1<dim; d1++) {
1156 for (
int d2=d1; d2<dim; d2++) {
1157 dealii::Point<dim,real> pert_p_p = point;
1158 dealii::Point<dim,real> pert_p_m = point;
1159 dealii::Point<dim,real> pert_m_p = point;
1160 dealii::Point<dim,real> pert_m_m = point;
1162 pert_p_p[d1] += (+eps); pert_p_p[d2] += (+eps);
1163 pert_p_m[d1] += (+eps); pert_p_m[d2] += (-eps);
1164 pert_m_p[d1] += (-eps); pert_m_p[d2] += (+eps);
1165 pert_m_m[d1] += (-eps); pert_m_m[d2] += (-eps);
1167 const real valpp = value(pert_p_p, istate);
1168 const real valpm = value(pert_p_m, istate);
1169 const real valmp = value(pert_m_p, istate);
1170 const real valmm = value(pert_m_m, istate);
1172 hessian[d1][d2] = (valpp - valpm - valmp + valmm) / (4*eps*eps);
1178 template <
int dim,
typename real>
1181 const dealii::Point<dim,real> &p,
1182 std::vector<dealii::Tensor<1,dim, real> > &gradients)
const 1184 for (
unsigned int i = 0; i < nstate; ++i)
1185 gradients[i] = gradient(p, i);
1189 template <
int dim,
typename real>
1193 std::vector<real> values(nstate);
1194 for (
unsigned int s=0; s<nstate; s++) { values[s] = value(point, s); }
1198 template <
int dim,
typename real>
1199 std::shared_ptr< ManufacturedSolutionFunction<dim,real> >
1207 return create_ManufacturedSolution(solution_type, nstate);
1210 template <
int dim,
typename real>
1211 std::shared_ptr< ManufacturedSolutionFunction<dim,real> >
1216 if(solution_type == ManufacturedSolutionEnum::sine_solution){
1217 return std::make_shared<ManufacturedSolutionSine<dim,real>>(nstate);
1218 }
else if(solution_type == ManufacturedSolutionEnum::zero_solution){
1219 return std::make_shared<ManufacturedSolutionZero<dim,real>>(nstate);
1220 }
else if(solution_type == ManufacturedSolutionEnum::cosine_solution){
1221 return std::make_shared<ManufacturedSolutionCosine<dim,real>>(nstate);
1222 }
else if(solution_type == ManufacturedSolutionEnum::additive_solution){
1223 return std::make_shared<ManufacturedSolutionAdd<dim,real>>(nstate);
1224 }
else if(solution_type == ManufacturedSolutionEnum::exp_solution){
1225 return std::make_shared<ManufacturedSolutionExp<dim,real>>(nstate);
1226 }
else if(solution_type == ManufacturedSolutionEnum::poly_solution){
1227 return std::make_shared<ManufacturedSolutionPoly<dim,real>>(nstate);
1228 }
else if(solution_type == ManufacturedSolutionEnum::even_poly_solution){
1229 return std::make_shared<ManufacturedSolutionEvenPoly<dim,real>>(nstate);
1230 }
else if(solution_type == ManufacturedSolutionEnum::atan_solution){
1231 return std::make_shared<ManufacturedSolutionAtan<dim,real>>(nstate);
1232 }
else if(solution_type == ManufacturedSolutionEnum::boundary_layer_solution){
1233 return std::make_shared<ManufacturedSolutionBoundaryLayer<dim,real>>(nstate);
1234 }
else if((solution_type == ManufacturedSolutionEnum::s_shock_solution) && (dim==2)){
1235 return std::make_shared<ManufacturedSolutionSShock<dim,real>>(nstate);
1236 }
else if(solution_type == ManufacturedSolutionEnum::quadratic_solution){
1237 return std::make_shared<ManufacturedSolutionQuadratic<dim,real>>(nstate);
1238 }
else if(solution_type == ManufacturedSolutionEnum::example_solution){
1239 return std::make_shared<ManufacturedSolutionExample<dim,real>>(nstate);
1240 }
else if((solution_type == ManufacturedSolutionEnum::navah_solution_1) && (dim==2) && (nstate==dim+2 || nstate==dim+3)){
1241 return std::make_shared<ManufacturedSolutionNavah_MS1<dim,real>>(nstate);
1242 }
else if((solution_type == ManufacturedSolutionEnum::navah_solution_2) && (dim==2) && (nstate==dim+2 || nstate==dim+3)){
1243 return std::make_shared<ManufacturedSolutionNavah_MS2<dim,real>>(nstate);
1244 }
else if((solution_type == ManufacturedSolutionEnum::navah_solution_3) && (dim==2) && (nstate==dim+2 || nstate==dim+3)){
1245 return std::make_shared<ManufacturedSolutionNavah_MS3<dim,real>>(nstate);
1246 }
else if((solution_type == ManufacturedSolutionEnum::navah_solution_4) && (dim==2) && (nstate==dim+2 || nstate==dim+3)){
1247 return std::make_shared<ManufacturedSolutionNavah_MS4<dim,real>>(nstate);
1248 }
else if((solution_type == ManufacturedSolutionEnum::navah_solution_5) && (dim==2) && (nstate==dim+2 || nstate==dim+3)){
1249 return std::make_shared<ManufacturedSolutionNavah_MS5<dim,real>>(nstate);
1251 std::cout <<
"Invalid combination of Manufactured Solution, dimension, and PDE Type." << std::endl;
1256 using FadType = Sacado::Fad::DFad<double>;
1257 using FadFadType = Sacado::Fad::DFad<FadType>;
1262 using codi_FadType = codi::RealForwardGen<double, codi::Direction<double,dimForwardAD>>;
1267 codi::Direction< codi::RealForwardVec<dimForwardAD>, dimReverseAD>
ManufacturedSolutionType
Selects the manufactured solution to be used if use_manufactured_source_term=true.
real value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Value.
Sum of even order polynomial functions manufactured solution.
codi::RealReverseIndexVec< dimReverseAD > codi_JacobianComputationType
Reverse mode type for Jacobian computation using TapeHelper.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Hessian.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Gradient.
Sacado::Fad::DFad< FadType > FadFadType
Sacado AD type that allows 2nd derivatives.
real value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Value.
Navah and Nadarajah free flows manufactured solution: MS5.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Hessian.
Sacado::Fad::DFad< double > FadType
Sacado AD type for first derivatives.
real value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Value.
Manufactured solution used for grid studies to check convergence orders.
codi_JacobianComputationType RadType
CoDiPaco reverse-AD type for first derivatives.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Hessian.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Gradient.
static std::shared_ptr< ManufacturedSolutionFunction< dim, real > > create_ManufacturedSolution(Parameters::AllParameters const *const param, int nstate)
Construct Manufactured solution object from global parameter file.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Gradient.
Files for the baseline physics.
Manufactured solution function factory.
real primitive_value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const
Value of primitive variables.
ManufacturedSolutionParam manufactured_solution_param
Associated manufactured solution parameters.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Hessian.
real value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Value.
real value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Value.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Hessian.
Main parameter class that contains the various other sub-parameter classes.
real value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Value.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Gradient.
Sum of polynomial manufactured solution.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &, const unsigned int=0) const override
Gradient.
ManufacturedConvergenceStudyParam manufactured_convergence_study_param
Contains parameters for manufactured convergence study.
Scalar boundary layer manufactured solution.
Product of cosine waves manufactured solution.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Gradient of conservative variables.
real value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Value of conservative variables.
static constexpr int dimForwardAD
Size of the forward vector mode for CoDiPack.
Navah and Nadarajah free flows manufactured solution: MS3.
std::vector< real > stdvector_values(const dealii::Point< dim, real > &point) const
Same as Function::values() except it returns it into a std::vector format.
Sum of exponential functions manufactured solution.
real value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Value.
dealii::SymmetricTensor< 2, dim, real > hessian_fd(const dealii::Point< dim, real > &point, const unsigned int istate=0) const
Uses finite-difference to evaluate the hessian.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Hessian of conservative variables.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Hessian.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Gradient.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &, const unsigned int=0) const override
Hessian.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Gradient.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Hessian.
Navah and Nadarajah free flows manufactured solution: MS1.
Product of sine waves manufactured solution.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Gradient.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Hessian.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Hessian.
codi::RealReversePrimalIndexGen< codi::RealForwardVec< dimForwardAD >, codi::Direction< codi::RealForwardVec< dimForwardAD >, dimReverseAD > > codi_HessianComputationType
Nested reverse-forward mode type for Jacobian and Hessian computation using TapeHelper.
Navah and Nadarajah free flows manufactured solution: MS4.
Hump manufactured solution based on arctangent functions.
real value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Value.
bool isfinite(double value)
< Provide isfinite for double.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Hessian.
Sum of sine waves manufactured solution.
dealii::SymmetricTensor< 2, dim, real > hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Hessian.
real value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Value.
dealii::Tensor< 1, dim, real > primitive_gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const
Gradient of primitive variables.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Gradient.
dealii::SymmetricTensor< 2, dim, real > primitive_hessian(const dealii::Point< dim, real > &point, const unsigned int istate=0) const
Hessian of primitive variables.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Gradient.
Product of zero waves manufactured solution.
ManufacturedSolutionType manufactured_solution_type
Selected ManufacturedSolutionType from the input file.
void vector_gradient(const dealii::Point< dim, real > &p, std::vector< dealii::Tensor< 1, dim, real > > &gradients) const
See dealii::Function<dim,real>::vector_gradient.
static constexpr int dimReverseAD
Size of the reverse vector mode for CoDiPack.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Gradient.
codi::RealForwardGen< double, codi::Direction< double, dimForwardAD > > codi_FadType
Tapeless forward mode.
ManufacturedSolutionFunction(const unsigned int nstate=1)
Constructor that initializes base_values, amplitudes, frequencies.
Navah and Nadarajah free flows manufactured solution: MS2.
S-Shock manufactured solution.
codi_HessianComputationType RadFadType
Nested reverse-forward mode type for Jacobian and Hessian computation using TapeHelper.
real value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Value.
dealii::Tensor< 1, dim, real > gradient_fd(const dealii::Point< dim, real > &point, const unsigned int istate=0) const
Uses finite-difference to evaluate the gradient.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Gradient.
Quadratic function manufactured solution.
real value(const dealii::Point< dim, real > &, const unsigned int=0) const override
Value.
real value(const dealii::Point< dim, real > &point, const unsigned int istate=0) const override
Value.