[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
manufactured_solution.cpp
1 #include <CoDiPack/include/codi.hpp>
2 #include <Sacado.hpp>
3 #include <deal.II/base/function.h>
4 #include <deal.II/base/function.templates.h> // Needed to instantiate dealii::Function<PHILIP_DIM,Sacado::Fad::DFad<double>>
5 #include <deal.II/base/function_time.templates.h> // Needed to instantiate dealii::Function<PHILIP_DIM,Sacado::Fad::DFad<double>>
6 
7 #include "manufactured_solution.h"
8 
9 template class dealii::FunctionTime<Sacado::Fad::DFad<double>>; // Needed by Function
10 template class dealii::Function<PHILIP_DIM,Sacado::Fad::DFad<double>>;
11 
12 namespace PHiLiP {
13 
15 bool isfinite(double value)
16 {
17  return std::isfinite(static_cast<double>(value));
18 }
19 
21 bool isfinite(Sacado::Fad::DFad<double> value)
22 {
23  return std::isfinite(static_cast<double>(value.val()));
24 }
25 
27 bool isfinite(Sacado::Fad::DFad<Sacado::Fad::DFad<double>> value)
28 {
29  return std::isfinite(static_cast<double>(value.val().val()));
30 }
31 
33 bool isfinite(Sacado::Rad::ADvar<Sacado::Fad::DFad<double>> value)
34 {
35  return std::isfinite(static_cast<double>(value.val().val()));
36 }
37 
38 template <int dim, typename real>
40 ::value (const dealii::Point<dim,real> &/*point*/, const unsigned int /*istate*/) const
41 {
42  real value = 0;
43  return value;
44 }
45 
46 template <int dim, typename real>
48 ::value (const dealii::Point<dim,real> &point, const unsigned int istate) const
49 {
50  real value = this->amplitudes[istate];
51  for (int d=0; d<dim; d++) {
52  value *= sin( this->frequencies[istate][d] * point[d] );
53  assert(isfinite(value));
54  }
55  value += this->base_values[istate];
56  return value;
57 }
58 
59 template <int dim, typename real>
61 ::value (const dealii::Point<dim,real> &point, const unsigned int istate) const
62 {
63  real value = 0.0;
64  for (int d=0; d<dim; d++) {
65  value += this->amplitudes[istate]*sin( this->frequencies[istate][d] * point[d] );
66  assert(isfinite(value));
67  }
68  value += this->base_values[istate];
69  return value;
70 }
71 
72 template <int dim, typename real>
74 ::value (const dealii::Point<dim,real> &point, const unsigned int istate) const
75 {
76  real value = this->amplitudes[istate];
77  for (int d=0; d<dim; d++) {
78  value *= cos( this->frequencies[istate][d] * point[d] );
79  assert(isfinite(value));
80  }
81  value += this->base_values[istate];
82  return value;
83 }
84 
85 template <int dim, typename real>
87 ::value (const dealii::Point<dim,real> &point, const unsigned int istate) const
88 {
89  real value = 0.0;
90  for (int d=0; d<dim; d++) {
91  value += exp( point[d] );
92  assert(isfinite(value));
93  }
94  value += this->base_values[istate];
95  return value;
96 }
97 
98 template <int dim, typename real>
100 ::value (const dealii::Point<dim,real> &point, const unsigned int istate) const
101 {
102  real value = 0.0;
103  const double poly_max = 7;
104  for (int d=0; d<dim; d++) {
105  value += pow(point[d] + 0.5, poly_max);
106  }
107  value += this->base_values[istate];
108  return value;
109 }
110 
111 template <int dim, typename real>
113 ::value (const dealii::Point<dim,real> &point, const unsigned int istate) const
114 {
115  real value = 0.0;
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);
119  }
120  value += this->base_values[istate];
121  return value;
122 }
123 
124 template <int dim, typename real>
126 ::value(const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
127 {
128  real val = 1.0;
129  for(unsigned int i = 0; i < dim; ++i){
130  real x = point[i];
131  real val_dim = 0;
132  for(unsigned int j = 0; j < n_shocks[i]; ++j){
133  // taking the product of function in each direction
134  val_dim += atan(S_j[i][j]*(x-x_j[i][j]));
135  }
136  val *= val_dim;
137  }
138  return val;
139 }
140 
141 template <int dim, typename real>
143 ::value(const dealii::Point<dim,real> &point, const unsigned int istate) const
144 {
145  real val = 1.0;
146  for(unsigned int d = 0; d < dim; ++d){
147  real x = point[d];
148  val *= x + (exp(x/epsilon[istate][d])-1.0)/(1.0-exp(1.0/epsilon[istate][d]));
149  }
150  return val;
151 }
152 
153 template <int dim, typename real>
155 ::value(const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
156 {
157  real val = 0.0;
158  if(dim==2){
159  const real x = point[0], y = point[1];
160  // val = 0.75*tanh(2*(sin(5*y)-3*x));
161  // val = 0.75*tanh(20*(sin(10*y-5)-6*x+3));
162  val = a*tanh(b*sin(c*y + d) + e*x + f);
163  }
164  return val;
165 }
166 
167 template <int dim, typename real>
169 ::value(const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
170 {
171  real val = 0.0;
172  // f(x,y,z) = a*x^2 + b*y^2 + c*z^2
173  for(unsigned int d = 0; d < dim; ++d){
174  real x = point[d];
175  val += alpha_diag[d]*x*x;
176  }
177  return val;
178 }
179 
180 template <int dim, typename real>
182 ::value (const dealii::Point<dim,real> &point, const unsigned int istate) const
183 {
184  real value = 0.0;
185  for (int d=0; d<dim; d++) {
186  value += exp( point[d] ) + sin(point[d] );
187  assert(isfinite(value));
188  }
189  value += this->base_values[istate];
190  return value;
191 }
192 
193 template <int dim, typename real>
195 ::primitive_value(const dealii::Point<dim,real> &point, const unsigned int istate) const
196 {
197  real value = 0.;
198  if constexpr(dim == 2) {
199  const real x = point[0], y = point[1];
200  // // for RANS
201  // const real v_tilde = ncm[istate][0] + ncm[istate][1]*cos(ncm[istate][4]*c*x) + ncm[istate][2]*cos(ncm[istate][5]*c*y) + ncm[istate][3]*cos(ncm[istate][6]*c*x)*cos(ncm[istate][6]*c*y);
202 
203  if(istate==0) {
204  // density
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);
206  }
207  if(istate==1) {
208  // x-velocity
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);
210  }
211  if(istate==2) {
212  // y-velocity
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);
214  }
215  if(istate==3) {
216  // pressure
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);
218  }
219  if(istate==4) {
220  // turbulent working variable
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);
222  }
223  }
224  return value;
225 }
226 
227 template <int dim, typename real>
229 ::value(const dealii::Point<dim,real> &point, const unsigned int istate) const
230 {
231  real value = 0.0;
232  if (dim == 2) {
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);
238 
239  // convert primitive to conservative solution
240  // - density:
241  if(istate==0) value = density;
242  // - x-momentum:
243  if(istate==1) value = density*x_velocity;
244  // - y-momentum:
245  if(istate==2) value = density*y_velocity;
246  // - total energy:
247  if(istate==3) value = pressure/(1.4-1.0) + 0.5*density*(x_velocity*x_velocity + y_velocity*y_velocity);
248  // - transport of the turbulent working variable:
249  if(istate==4) value = density*turbulent_working_variable;
250  }
251  return value;
252 }
253 
254 template <int dim, typename real>
255 inline dealii::Tensor<1,dim,real> ManufacturedSolutionZero<dim,real>
256 ::gradient (const dealii::Point<dim,real> &/*point*/, const unsigned int /*istate*/) const
257 {
258  dealii::Tensor<1,dim,real> gradient;
259  for(unsigned int i = 0; i < dim; i++){
260  gradient[i] = 0;
261  }
262  return gradient;
263 }
264 
265 template <int dim, typename real>
266 inline dealii::Tensor<1,dim,real> ManufacturedSolutionSine<dim,real>
267 ::gradient (const dealii::Point<dim,real> &point, const unsigned int istate) const
268 {
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 );
276  }
277  assert(isfinite(gradient[dim_deri]));
278  }
279  // Hard-coded is much more readable than the dimensionally generic one
280  const real A = this->amplitudes[istate];
281  const dealii::Tensor<1,dim,real> f = this->frequencies[istate];
282  if (dim==1) {
283  const real fx = f[0]*point[0];
284  gradient[0] = A*f[0]*cos(fx);
285  }
286  if (dim==2) {
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);
291  }
292  if (dim==3) {
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);
299  }
300  return gradient;
301 }
302 
303 template <int dim, typename real>
304 inline dealii::Tensor<1,dim,real> ManufacturedSolutionAdd<dim,real>
305 ::gradient (const dealii::Point<dim,real> &point, const unsigned int istate) const
306 {
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];
310  if (dim==1) {
311  const real fx = f[0]*point[0];
312  gradient[0] = A*f[0]*cos(fx);
313  }
314  if (dim==2) {
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);
319  }
320  if (dim==3) {
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);
327  }
328  return gradient;
329 }
330 
331 template <int dim, typename real>
332 inline dealii::Tensor<1,dim,real> ManufacturedSolutionCosine<dim,real>
333 ::gradient (const dealii::Point<dim,real> &point, const unsigned int istate) const
334 {
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];
338  if (dim==1) {
339  const real fx = f[0]*point[0];
340  gradient[0] = -A*f[0]*sin(fx);
341  }
342  if (dim==2) {
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);
347  }
348  if (dim==3) {
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);
355  }
356  return gradient;
357 }
358 
359 template <int dim, typename real>
360 inline dealii::Tensor<1,dim,real> ManufacturedSolutionExp<dim,real>
361 ::gradient (const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
362 {
363  dealii::Tensor<1,dim,real> gradient;
364  if (dim==1) {
365  gradient[0] = exp(point[0]);
366  }
367  if (dim==2) {
368  gradient[0] = exp(point[0]);
369  gradient[1] = exp(point[1]);
370  }
371  if (dim==3) {
372  gradient[0] = exp(point[0]);
373  gradient[1] = exp(point[1]);
374  gradient[2] = exp(point[2]);
375  }
376  return gradient;
377 }
378 
379 template <int dim, typename real>
380 inline dealii::Tensor<1,dim,real> ManufacturedSolutionEvenPoly<dim,real>
381 ::gradient (const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
382 {
383  dealii::Tensor<1,dim,real> gradient;
384  const double poly_max = 7;
385  if (dim==1) {
386  gradient[0] = poly_max*pow(point[0] + 0.5, poly_max-1);
387  }
388  if (dim==2) {
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);
391  }
392  if (dim==3) {
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);
396  }
397  return gradient;
398 }
399 
400 template <int dim, typename real>
401 inline dealii::Tensor<1,dim,real> ManufacturedSolutionPoly<dim,real>
402 ::gradient (const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
403 {
404  dealii::Tensor<1,dim,real> gradient;
405  if (dim==1) {
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);
408  }
409  if (dim==2) {
410  real x = point[0];
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);
412  x = point[1];
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);
414  }
415  if (dim==3) {
416  real x = point[0];
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;
418  x = point[1];
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;
420  x = point[2];
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;
422  }
423  return gradient;
424 }
425 
426 template <int dim, typename real>
427 inline dealii::Tensor<1,dim,real> ManufacturedSolutionAtan<dim,real>
428 ::gradient(const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
429 {
430  dealii::Tensor<1,dim,real> gradient;
431  for(unsigned int k = 0; k < dim; ++k){
432  // taking the k^th derivative
433  real grad_dim = 1;
434  for(unsigned int i = 0; i < dim; ++i){
435  real x = point[i];
436  real val_dim = 0;
437  for(unsigned int j = 0; j < n_shocks[i]; ++j){
438  if(i==k){
439  // taking the derivative dimension
440  real coeff = S_j[i][j]*(x-x_j[i][j]);
441  val_dim += S_j[i][j]/(pow(coeff,2)+1);
442  }else{
443  // value product unaffected
444  val_dim += atan(S_j[i][j]*(x-x_j[i][j]));
445  }
446  }
447  grad_dim *= val_dim;
448  }
449  gradient[k] = grad_dim;
450  }
451  return gradient;
452 }
453 
454 template <int dim, typename real>
455 inline dealii::Tensor<1,dim,real> ManufacturedSolutionBoundaryLayer<dim,real>
456 ::gradient(const dealii::Point<dim,real> &point, const unsigned int istate) const
457 {
458  dealii::Tensor<1,dim,real> gradient;
459  if(dim == 1){
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])));
462  }else if(dim == 2){
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])));
468  }else if(dim == 3){
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])));
479  }
480  return gradient;
481 }
482 
483 template <int dim, typename real>
484 inline dealii::Tensor<1,dim,real> ManufacturedSolutionSShock<dim,real>
485 ::gradient(const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
486 {
487  dealii::Tensor<1,dim,real> gradient;
488  if(dim == 2){
489  const real x = point[0], y = point[1];
490  // gradient[0] = -4.5*pow(cosh(6*x-2*sin(5*y)),-2);
491  // gradient[1] = 7.5*pow(cosh(6*x-2*sin(5*y)),-2)*cos(5*y);
492  // gradient[0] = -90*pow(cosh(-120*x-20*sin(5-10*y)+60),-2);
493  // gradient[1] = 150*pow(cosh(-120*x-20*sin(5-10*y)+60),-2)*cos(5-10*y);
494 
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;
498  }
499  return gradient;
500 }
501 
502 template <int dim, typename real>
503 inline dealii::Tensor<1,dim,real> ManufacturedSolutionQuadratic<dim,real>
504 ::gradient(const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
505 {
506  dealii::Tensor<1,dim,real> gradient;
507  for(unsigned int d = 0; d < dim; ++d){
508  // dF = <2ax, 2by, 2cz>
509  const real x = point[d];
510  gradient[d] = 2*alpha_diag[d]*x;
511  }
512  return gradient;
513 }
514 
515 template <int dim, typename real>
516 inline dealii::Tensor<1,dim,real> ManufacturedSolutionExample<dim,real>
517 ::gradient(const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
518 {
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]);
522  }
523  return gradient;
524 }
525 
526 template <int dim, typename real>
527 inline dealii::Tensor<1,dim,real> ManufacturedSolutionNavahBase<dim,real>
528 ::primitive_gradient (const dealii::Point<dim,real> &point, const unsigned int istate) const
529 {
530  dealii::Tensor<1,dim,real> gradient;
531  // Gradients of primitive variables
532  if (dim == 2) {
533  const real x = point[0], y = point[1];
534 
535  if(istate==0) {
536  // density
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); // dx
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); // dy
539  }
540  if(istate==1) {
541  // x-velocity
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); // dx
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); // dy
544  }
545  if(istate==2) {
546  // y-velocity
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); // dx
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); // dy
549  }
550  if(istate==3) {
551  // pressure
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); // dx
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); // dy
554  }
555  if(istate==4) {
556  // turbulent working variable
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); // dx
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); // dy
559  }
560  }
561  return gradient;
562 }
563 
564 template <int dim, typename real>
565 inline dealii::Tensor<1,dim,real> ManufacturedSolutionNavahBase<dim,real>
566 ::gradient (const dealii::Point<dim,real> &point, const unsigned int istate) const
567 {
568  dealii::Tensor<1,dim,real> gradient;
569 
570  if (dim == 2) {
571  const real rho = primitive_value(point,0);
572  const real u = primitive_value(point,1);
573  const real v = primitive_value(point,2);
574  // const real p = primitive_value(point,3);
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);
579 
580  // convert to primitive to gradient of conservative variables using product rule
581  if(istate==0) {
582  // density
583  for(int d=0; d<dim; d++) {
584  gradient[d] = rho_grad[d];
585  }
586  }
587  if(istate==1) {
588  // x-momentum
589  for(int d=0; d<dim; d++) {
590  gradient[d] = u*rho_grad[d] + rho*u_grad[d];
591  }
592  }
593  if(istate==2) {
594  // y-momentum
595  for(int d=0; d<dim; d++) {
596  gradient[d] = v*rho_grad[d] + rho*v_grad[d];
597  }
598  }
599  if(istate==3) {
600  // total energy
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]);
603  }
604  }
605  if(istate==4) {
606  // transport of the turbulent working variable (twv)
607  const real twv = primitive_value(point,4);
608  const dealii::Tensor<1,dim,real> twv_grad = primitive_gradient(point,4); // only used for RANS
609  for(int d=0; d<dim; d++) {
610  gradient[d] = twv*rho_grad[d] + rho*twv_grad[d];
611  }
612  }
613  }
614  return gradient;
615 }
616 
617 template <int dim, typename real>
618 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionZero<dim,real>
619 ::hessian (const dealii::Point<dim,real> &/*point*/, const unsigned int /*istate*/) const
620 {
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++){
624  hessian[i][j] = 0;
625  }
626  }
627  return hessian;
628 }
629 
630 template <int dim, typename real>
631 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionSine<dim,real>
632 ::hessian (const dealii::Point<dim,real> &point, const unsigned int istate) const
633 {
634  dealii::SymmetricTensor<2,dim,real> hessian;
635  // Hard-coded is much more readable than the dimensionally generic one
636  const real A = this->amplitudes[istate];
637  const dealii::Tensor<1,dim,real> f = this->frequencies[istate];
638  if (dim==1) {
639  const real fx = f[0]*point[0];
640  hessian[0][0] = -A*f[0]*f[0]*sin(fx);
641  }
642  if (dim==2) {
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);
647 
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);
650  }
651  if (dim==3) {
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);
658 
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);
662 
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);
666  }
667  return hessian;
668 }
669 
670 template <int dim, typename real>
671 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionAdd<dim,real>
672 ::hessian (const dealii::Point<dim,real> &point, const unsigned int istate) const
673 {
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];
677  if (dim==1) {
678  const real fx = f[0]*point[0];
679  hessian[0][0] = -A*f[0]*f[0]*sin(fx);
680  }
681  if (dim==2) {
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);
685  hessian[0][1] = 0.0;
686 
687  hessian[1][0] = 0.0;
688  hessian[1][1] = -A*f[1]*f[1]*sin(fy);
689  }
690  if (dim==3) {
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);
695  hessian[0][1] = 0.0;
696  hessian[0][2] = 0.0;
697 
698  hessian[1][0] = 0.0;
699  hessian[1][1] = -A*f[1]*f[1]*sin(fy);
700  hessian[1][2] = 0.0;
701 
702  hessian[2][0] = 0.0;
703  hessian[2][1] = 0.0;
704  hessian[2][2] = -A*f[2]*f[2]*sin(fz);
705  }
706  return hessian;
707 }
708 
709 template <int dim, typename real>
710 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionCosine<dim,real>
711 ::hessian (const dealii::Point<dim,real> &point, const unsigned int istate) const
712 {
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];
716  if (dim==1) {
717  const real fx = f[0]*point[0];
718  hessian[0][0] = -A*f[0]*f[0]*cos(fx);
719  }
720  if (dim==2) {
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);
725 
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);
728  }
729  if (dim==3) {
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);
736 
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);
740 
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);
744  }
745  return hessian;
746 }
747 
748 template <int dim, typename real>
749 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionExp<dim,real>
750 ::hessian (const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
751 {
752  dealii::SymmetricTensor<2,dim,real> hessian;
753  if (dim==1) {
754  hessian[0][0] = exp(point[0]);
755  }
756  if (dim==2) {
757  hessian[0][0] = exp(point[0]);
758  hessian[0][1] = 0.0;
759 
760  hessian[1][0] = 0.0;
761  hessian[1][1] = exp(point[1]);
762  }
763  if (dim==3) {
764  hessian[0][0] = exp(point[0]);
765  hessian[0][1] = 0.0;
766  hessian[0][2] = 0.0;
767 
768  hessian[1][0] = 0.0;
769  hessian[1][1] = exp(point[1]);
770  hessian[1][2] = 0.0;
771 
772  hessian[2][0] = 0.0;
773  hessian[2][1] = 0.0;
774  hessian[2][2] = exp(point[2]);
775  }
776  return hessian;
777 }
778 
779 template <int dim, typename real>
780 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionEvenPoly<dim,real>
781 ::hessian (const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
782 {
783  dealii::SymmetricTensor<2,dim,real> hessian;
784  const double poly_max = 7;
785  if (dim==1) {
786  hessian[0][0] = poly_max*poly_max*pow(point[0] + 0.5, poly_max-2);
787  }
788  if (dim==2) {
789  hessian[0][0] = poly_max*poly_max*pow(point[0] + 0.5, poly_max-2);
790  hessian[0][1] = 0.0;
791 
792  hessian[1][0] = 0.0;
793  hessian[1][1] = poly_max*poly_max*pow(point[1] + 0.5, poly_max-2);
794  }
795  if (dim==3) {
796  hessian[0][0] = poly_max*poly_max*pow(point[0] + 0.5, poly_max-2);
797  hessian[0][1] = 0.0;
798  hessian[0][2] = 0.0;
799 
800  hessian[1][0] = 0.0;
801  hessian[1][1] = poly_max*poly_max*pow(point[1] + 0.5, poly_max-2);
802  hessian[1][2] = 0.0;
803 
804  hessian[2][0] = 0.0;
805  hessian[2][1] = 0.0;
806  hessian[2][2] = poly_max*poly_max*pow(point[2] + 0.5, poly_max-2);
807  }
808  return hessian;
809 }
810 
811 template <int dim, typename real>
812 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionPoly<dim,real>
813 ::hessian (const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
814 {
815  dealii::SymmetricTensor<2,dim,real> hessian;
816  if (dim==1) {
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);
819  }
820  if (dim==2) {
821  real x = point[0];
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);
823  x = point[1];
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);
825  }
826  if (dim==3) {
827  real x = point[0];
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);
829  x = point[1];
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);
831  x = point[2];
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);
833  }
834  return hessian;
835 }
836 
837 template <int dim, typename real>
838 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionAtan<dim,real>
839 ::hessian(const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
840 {
841  dealii::SymmetricTensor<2,dim,real> hes;
842 
843  for(unsigned int k1 = 0; k1 < dim; ++k1){
844  // taking the k1^th derivative
845  for(unsigned int k2 = 0; k2 < dim; ++k2){
846  // taking the k2^th derivative
847  real hes_dim = 1;
848  for(unsigned int i = 0; i < dim; ++i){
849  real x = point[i];
850  real val_dim = 0;
851  for(unsigned int j = 0; j < n_shocks[i]; ++j){
852  if(i == k1 && i == k2){
853  // taking the second derivative in this dim
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){
857  // taking the first derivative in this dim
858  real coeff = S_j[i][j]*(x-x_j[i][j]);
859  val_dim += S_j[i][j]/(pow(coeff,2)+1);
860  }else{
861  // taking the value in this dim
862  val_dim += atan(S_j[i][j]*(x-x_j[i][j]));
863  }
864  }
865  hes_dim *= val_dim;
866  }
867  hes[k1][k2] = hes_dim;
868  }
869  }
870 
871  return hes;
872 }
873 
874 template <int dim, typename real>
875 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionBoundaryLayer<dim,real>
876 ::hessian(const dealii::Point<dim,real> &point, const unsigned int istate) const
877 {
878  dealii::SymmetricTensor<2,dim,real> hessian;
879  if (dim==1) {
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])));
882  }
883  if (dim==2) {
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])));
889 
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])));
893  }
894  if (dim==3) {
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])));
905 
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])));
913 
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])));
919  }
920  return hessian;
921 }
922 
923 template <int dim, typename real>
924 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionSShock<dim,real>
925 ::hessian(const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
926 {
927  dealii::SymmetricTensor<2,dim,real> hessian;
928  if (dim==2) {
929  const real x = point[0], y = point[1];
930  // hessian[0][0] = 54*tanh(6*x-2*sin(5*y))*pow(cosh(6*x-2*sin(5*y)),-2);
931  // hessian[0][1] = -90*tanh(6*x-2*sin(5*y))*pow(cosh(6*x-2*sin(5*y)),-2)*cos(5*y);
932 
933  // hessian[1][0] = hessian[1][0];
934  // hessian[1][1] = pow(cosh(6*x-2*sin(5*y)),-2)*(-37.5*sin(5*y)+150*pow(cos(5*y),2)*tanh(6*x-2*sin(5*y)));
935 
936  // hessian[0][0] = 21600*pow(cosh(20*(-3+6*x+sin(5-10*y))),2)*tanh(20*(-3+6*x+sin(5-10*y)));
937  // hessian[0][1] = -36000*pow(cosh(20*(-3+6*x+sin(5-10*y))),2)*tanh(20*(-3+6*x+sin(5-10*y)))*cos(5-10*y);
938 
939  // hessian[1][0] = hessian[0][1];
940  // hessian[1][1] = 1500*pow(cosh(20*(-3+6*x+sin(5-10*y))),2)*(40*pow(cos(5-10*y),2)*tanh(20*(-3+6*x+sin(5-10*y)))+sin(5-10*y));
941 
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);
945 
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;
948 
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));
951  }
952  return hessian;
953 }
954 
955 template <int dim, typename real>
956 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionQuadratic<dim,real>
957 ::hessian(const dealii::Point<dim,real> &/* point */, const unsigned int /* istate */) const
958 {
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){
962  if(i == j){
963  hessian[i][i] = 2*alpha_diag[i];
964  }else{
965  hessian[i][j] = 0.0;
966  }
967  }
968  }
969  return hessian;
970 }
971 
972 template <int dim, typename real>
973 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionExample<dim,real>
974 ::hessian (const dealii::Point<dim,real> &point, const unsigned int /*istate*/) const
975 {
976  dealii::SymmetricTensor<2,dim,real> hessian;
977  for(int idim=0; idim<dim; idim++){
978  for(int jdim=0; jdim<dim; jdim++){
979  if(idim == jdim)
980  hessian[idim][jdim] = exp(point[idim]) - sin(point[idim]);
981  else
982  hessian[idim][jdim] = 0.0;
983  }
984  }
985  return hessian;
986 }
987 
988 template <int dim, typename real>
989 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionNavahBase<dim,real>
990 ::primitive_hessian (const dealii::Point<dim,real> &point, const unsigned int istate) const
991 {
992  dealii::SymmetricTensor<2,dim,real> hessian;
993 
994  if (dim == 2) {
995  const real x = point[0], y = point[1];
996 
997  if(istate==0) {
998  // density
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); // dxdx
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); // dxdy
1001  hessian[1][0] = hessian[0][1]; // dydx
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); // dydy
1003  }
1004  if(istate==1) {
1005  // x-velocity
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); // dxdx
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); // dxdy
1008  hessian[1][0] = hessian[0][1]; // dydx
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); // dydy
1010  }
1011  if(istate==2) {
1012  // y-velocity
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); // dxdx
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); // dxdy
1015  hessian[1][0] = hessian[0][1]; // dydx
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); // dydy
1017  }
1018  if(istate==3) {
1019  // pressure
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); // dxdx
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); // dxdy
1022  hessian[1][0] = hessian[0][1]; // dydx
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); // dydy
1024  }
1025  if(istate==4) {
1026  // turbulent working variable
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); // dxdx
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); // dxdy
1029  hessian[1][0] = hessian[0][1]; // dydx
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); // dydy
1031  }
1032  }
1033  return hessian;
1034 }
1035 
1036 template <int dim, typename real>
1037 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionNavahBase<dim,real>
1038 ::hessian (const dealii::Point<dim,real> &point, const unsigned int istate) const
1039 {
1040  dealii::SymmetricTensor<2,dim,real> hessian;
1041 
1042  if (dim == 2) {
1043  const real rho = primitive_value(point,0);
1044  const real u = primitive_value(point,1);
1045  const real v = primitive_value(point,2);
1046  // const real p = primitive_value(point,3);
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);
1050  // const dealii::Tensor<1,dim,real> p_grad = primitive_gradient(point,3);
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);
1055 
1056  // convert to primitive to hessian of conservative variables using product rule
1057  if(istate==0) {
1058  // density
1059  for(int i=0; i<dim; i++) {
1060  for(int j=0; j<dim; j++) {
1061  hessian[i][j] = rho_hess[i][j];
1062  }
1063  }
1064  }
1065  if(istate==1) {
1066  // x-momentum
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];
1070  }
1071  }
1072  }
1073  if(istate==2) {
1074  // y-momentum
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];
1078  }
1079  }
1080  }
1081  if(istate==3) {
1082  // total energy
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]);
1088  }
1089  }
1090  }
1091  if(istate==4) {
1092  // transport of the turbulent working variable (twv)
1093  const real twv = primitive_value(point,4); // only used by istate=4
1094  const dealii::Tensor<1,dim,real> twv_grad = primitive_gradient(point,4); // only used by istate=4
1095  const dealii::SymmetricTensor<2,dim,real> twv_hess = primitive_hessian(point,4); // only used by istate=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];
1099  }
1100  }
1101  }
1102  }
1103  return hessian;
1104 }
1105 
1106 template <int dim, typename real>
1108 ::ManufacturedSolutionFunction (const unsigned int nstate)
1109  :
1110  dealii::Function<dim,real>(nstate)
1111  , nstate(nstate)
1112  , base_values(nstate)
1113  , amplitudes(nstate)
1114  , frequencies(nstate)
1115 {
1116  const double pi = atan(1)*4.0;
1117  //const double ee = exp(1);
1118 
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++) {
1124  //frequencies[s][d] = 2.0 + sin(0.1+s*0.5+d*0.2) * pi / 2.0;
1125  frequencies[s][d] = 2.0 + sin(0.1+s*0.5+d*0.2) * pi / 2.0;
1126  }
1127 
1128  }
1129 }
1130 
1131 template <int dim, typename real>
1132 inline dealii::Tensor<1,dim,real> ManufacturedSolutionFunction<dim,real>
1133 ::gradient_fd (const dealii::Point<dim,real> &point, const unsigned int istate) const
1134 {
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);
1145  }
1146  return gradient;
1147 }
1148 
1149 template <int dim, typename real>
1150 inline dealii::SymmetricTensor<2,dim,real> ManufacturedSolutionFunction<dim,real>
1151 ::hessian_fd (const dealii::Point<dim,real> &point, const unsigned int istate) const
1152 {
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;
1161 
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);
1166 
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);
1171 
1172  hessian[d1][d2] = (valpp - valpm - valmp + valmm) / (4*eps*eps);
1173  }
1174  }
1175  return hessian;
1176 }
1177 
1178 template <int dim, typename real>
1181  const dealii::Point<dim,real> &p,
1182  std::vector<dealii::Tensor<1,dim, real> > &gradients) const
1183 {
1184  for (unsigned int i = 0; i < nstate; ++i)
1185  gradients[i] = gradient(p, i);
1186 }
1187 
1188 
1189 template <int dim, typename real>
1190 inline std::vector<real> ManufacturedSolutionFunction<dim,real>
1191 ::stdvector_values (const dealii::Point<dim,real> &point) const
1192 {
1193  std::vector<real> values(nstate);
1194  for (unsigned int s=0; s<nstate; s++) { values[s] = value(point, s); }
1195  return values;
1196 }
1197 
1198 template <int dim, typename real>
1199 std::shared_ptr< ManufacturedSolutionFunction<dim,real> >
1201  Parameters::AllParameters const *const param,
1202  int nstate)
1203 {
1206 
1207  return create_ManufacturedSolution(solution_type, nstate);
1208 }
1209 
1210 template <int dim, typename real>
1211 std::shared_ptr< ManufacturedSolutionFunction<dim,real> >
1214  int nstate)
1215 {
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);
1250  }else{
1251  std::cout << "Invalid combination of Manufactured Solution, dimension, and PDE Type." << std::endl;
1252  }
1253  return nullptr;
1254 }
1255 
1256 using FadType = Sacado::Fad::DFad<double>;
1257 using FadFadType = Sacado::Fad::DFad<FadType>;
1258 
1259 static constexpr int dimForwardAD = 1;
1260 static constexpr int dimReverseAD = 1;
1261 
1262 using codi_FadType = codi::RealForwardGen<double, codi::Direction<double,dimForwardAD>>;
1263 //using codi_FadType = codi::RealForwardGen<double, codi::DirectionVar<double>>;
1264 
1265 using codi_JacobianComputationType = codi::RealReverseIndexVec<dimReverseAD>;
1266 using codi_HessianComputationType = codi::RealReversePrimalIndexGen< codi::RealForwardVec<dimForwardAD>,
1267  codi::Direction< codi::RealForwardVec<dimForwardAD>, dimReverseAD>
1268  >;
1269 //using RadFadType = Sacado::Rad::ADvar<FadType>; ///< Sacado AD type that allows 2nd derivatives.
1270 //using RadFadType = codi_JacobianComputationType; ///< Reverse only mode that only allows Jacobian computation.
1273 
1279 
1365 
1371 }
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.
Definition: ADTypes.hpp:20
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.
Definition: ADTypes.hpp:12
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.
Definition: ADTypes.hpp:11
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.
Definition: ADTypes.hpp:27
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.
Definition: ADTypes.hpp:10
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.
Definition: ADTypes.hpp:14
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.
Definition: ADTypes.hpp:23
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.
Definition: ADTypes.hpp:15
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.
Definition: ADTypes.hpp:17
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.
Definition: ADTypes.hpp:28
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.