[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
euler_cylinder_adjoint.cpp
1 #include "euler_cylinder_adjoint.h"
2 
3 #include <deal.II/base/convergence_table.h>
4 #include <deal.II/distributed/grid_refinement.h>
5 #include <deal.II/distributed/solution_transfer.h>
6 #include <deal.II/distributed/tria.h>
7 #include <deal.II/dofs/dof_tools.h>
8 #include <deal.II/fe/fe_values.h>
9 #include <deal.II/fe/mapping_q.h>
10 #include <deal.II/grid/grid_generator.h>
11 #include <deal.II/grid/grid_in.h>
12 #include <deal.II/grid/grid_out.h>
13 #include <deal.II/grid/grid_refinement.h>
14 #include <deal.II/grid/grid_tools.h>
15 #include <deal.II/numerics/vector_tools.h>
16 #include <stdlib.h>
17 
18 #include <iostream>
19 
20 #include "dg/dg_base.hpp"
21 #include "dg/dg_factory.hpp"
22 #include "functional/adjoint.h"
23 #include "functional/functional.h"
24 #include "ode_solver/ode_solver_factory.h"
25 #include "physics/euler.h"
26 #include "physics/initial_conditions/initial_condition_function.h"
27 #include "physics/manufactured_solution.h"
28 #include "physics/physics.h"
29 #include "physics/physics_factory.h"
30 
31 namespace PHiLiP {
32 namespace Tests {
33 
36 template <int dim, int nstate, typename real>
37 class L2normError : public Functional<dim, nstate, real>
38 {
39 public:
42  std::shared_ptr<PHiLiP::DGBase<dim,real>> dg_input,
43  const bool uses_solution_values = true,
44  const bool uses_solution_gradient = true)
45  : PHiLiP::Functional<dim,nstate,real>(dg_input,uses_solution_values,uses_solution_gradient)
46  {}
47 
49  template <typename real2>
52  const dealii::Point<dim,real2> &/*phys_coord*/,
53  const std::array<real2,nstate> &soln_at_q,
54  const std::array<dealii::Tensor<1,dim,real2>,nstate> &/*soln_grad_at_q*/) const
55  {
56 
57  real2 cell_l2error = 0;
58 
59  const Physics::Euler<dim,nstate,real2>& euler_physics = dynamic_cast<const Physics::Euler<dim,nstate,real2>&>(physics);
60 
61  const real2 entropy = euler_physics.compute_entropy_measure(soln_at_q);
62 
63  const real2 uexact = euler_physics.entropy_inf;
64  cell_l2error = pow(entropy - uexact, 2);
65 
66  return cell_l2error;
67  }
68 
72  const dealii::Point<dim,real> &phys_coord,
73  const std::array<real,nstate> &soln_at_q,
74  const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_at_q) const override
75  {
76  return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q);
77  }
78 
79  using FadType = Sacado::Fad::DFad<real>;
80  using FadFadType = Sacado::Fad::DFad<FadType>;
81 
85  const dealii::Point<dim,FadFadType> &phys_coord,
86  const std::array<FadFadType,nstate> &soln_at_q,
87  const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &soln_grad_at_q) const override
88  {
89  return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q);
90  }
91 };
92 
93 dealii::Point<2> center_adjoint(0.0,0.0);
94 const double inner_radius_adjoint = 1, outer_radius_adjoint = inner_radius_adjoint*20;
95 
96 dealii::Point<2> warp_cylinder_adjoint (const dealii::Point<2> &p)//, const double inner_radius_adjoint, const double outer_radius_adjoint)
97 {
98  const double rectangle_height = 1.0;
99  //const double original_radius = std::abs(p[0]);
100  const double angle = p[1]/rectangle_height * dealii::numbers::PI;
101 
102  //const double radius = std::abs(p[0]);
103 
104  const double power = 1.8;
105  const double radius = outer_radius_adjoint*(inner_radius_adjoint/outer_radius_adjoint + pow(std::abs(p[0]), power));
106 
107  dealii::Point<2> q = p;
108  q[0] = -radius*cos(angle);
109  q[1] = radius*sin(angle);
110  return q;
111 }
112 
113 void half_cylinder_adjoint(dealii::parallel::distributed::Triangulation<2> & tria,
114  const unsigned int n_cells_circle,
115  const unsigned int n_cells_radial)
116 {
117  //const double pi = dealii::numbers::PI;
118  //double inner_circumference = inner_radius_adjoint*pi;
119  //double outer_circumference = outer_radius_adjoint*pi;
120  //const double rectangle_height = inner_circumference;
121  dealii::Point<2> p1(-1,0.0), p2(-0.0,1.0);
122 
123  const bool colorize = true;
124 
125  std::vector<unsigned int> n_subdivisions(2);
126  n_subdivisions[0] = n_cells_radial;
127  n_subdivisions[1] = n_cells_circle;
128  dealii::GridGenerator::subdivided_hyper_rectangle (tria, n_subdivisions, p1, p2, colorize);
129 
130  dealii::GridTools::transform (&warp_cylinder_adjoint, tria);
131 
132  tria.set_all_manifold_ids(0);
133  tria.set_manifold(0, dealii::SphericalManifold<2>(center_adjoint));
134 
135  // Assign BC
136  for (auto cell = tria.begin_active(); cell != tria.end(); ++cell) {
137  //if (!cell->is_locally_owned()) continue;
138  for (unsigned int face=0; face<dealii::GeometryInfo<2>::faces_per_cell; ++face) {
139  if (cell->face(face)->at_boundary()) {
140  unsigned int current_id = cell->face(face)->boundary_id();
141  if (current_id == 0) {
142  cell->face(face)->set_boundary_id (1004); // x_left, Farfield
143  } else if (current_id == 1) {
144  cell->face(face)->set_boundary_id (1001); // x_right, Symmetry/Wall
145  } else if (current_id == 2) {
146  cell->face(face)->set_boundary_id (1001); // y_bottom, Symmetry/Wall
147  } else if (current_id == 3) {
148  cell->face(face)->set_boundary_id (1001); // y_top, Wall
149  } else {
150  std::abort();
151  }
152  }
153  }
154  }
155 }
156 
157 template <int dim, int nstate>
159  :
160  TestsBase::TestsBase(parameters_input)
161 {}
162 
163 template<int dim, int nstate>
165 ::run_test () const
166 {
167  pcout << " Running Euler cylinder adjoint entropy convergence. " << std::endl;
168  static_assert(dim==2);
170  using GridEnum = ManParam::GridEnum;
172 
173  Assert(dim == param.dimension, dealii::ExcDimensionMismatch(dim, param.dimension));
174  Assert(param.pde_type == param.PartialDifferentialEquation::euler, dealii::ExcMessage("Can't run Euler case if PDE is not Euler"));
175  //if (param.pde_type == param.PartialDifferentialEquation::euler) return 1;
176 
177  ManParam manu_grid_conv_param = param.manufactured_convergence_study_param;
178 
179  const unsigned int p_start = manu_grid_conv_param.degree_start;
180  const unsigned int p_end = manu_grid_conv_param.degree_end;
181 
182  const unsigned int n_grids_input = manu_grid_conv_param.number_of_grids;
183 
184  std::shared_ptr< Physics::PhysicsBase<dim,nstate,double> > physics_double
186 
187  std::shared_ptr< Physics::Euler<dim,nstate,double> > euler_physics_double
188  = std::dynamic_pointer_cast< Physics::Euler<dim,nstate,double> >(physics_double);
189 
190  std::shared_ptr< Physics::PhysicsBase<dim,nstate,Sacado::Fad::DFad<double>> > euler_physics_adtype
192 
193  // Physics::Euler<dim,nstate,double> euler_physics_double
194  // = Physics::Euler<dim, nstate, double>(
195  // param.euler_param.ref_length,
196  // param.euler_param.gamma_gas,
197  // param.euler_param.mach_inf,
198  // param.euler_param.angle_of_attack,
199  // param.euler_param.side_slip_angle);
200 
201  // Physics::Euler<dim,nstate,Sacado::Fad::DFad<double>> euler_physics_adtype
202  // = Physics::Euler<dim, nstate, Sacado::Fad::DFad<double>>(
203  // param.euler_param.ref_length,
204  // param.euler_param.gamma_gas,
205  // param.euler_param.mach_inf,
206  // param.euler_param.angle_of_attack,
207  // param.euler_param.side_slip_angle);
208 
209  FreeStreamInitialConditions<dim,nstate,double> initial_conditions(*euler_physics_double);
210 
211  std::vector<int> fail_conv_poly;
212  std::vector<double> fail_conv_slop;
213  std::vector<dealii::ConvergenceTable> convergence_table_vector;
214 
215  for (unsigned int poly_degree = p_start; poly_degree <= p_end; ++poly_degree) {
216 
217  // p0 tends to require a finer grid to reach asymptotic region
218  unsigned int n_grids = n_grids_input;
219  if (poly_degree <= 1) n_grids = n_grids_input;
220 
221  std::vector<double> entropy_error(n_grids);
222  std::vector<double> grid_size(n_grids);
223 
224  const std::vector<int> n_1d_cells = get_number_1d_cells(n_grids);
225 
226  dealii::ConvergenceTable convergence_table;
227 
228  // Generate grid and mapping
229  using Triangulation = dealii::parallel::distributed::Triangulation<dim>;
230  std::shared_ptr <Triangulation> grid = std::make_shared<Triangulation> (this->mpi_communicator);
231 
232  const unsigned int n_cells_circle = n_1d_cells[0];
233  const unsigned int n_cells_radial = 1.5*n_cells_circle;
234  half_cylinder_adjoint(*grid, n_cells_circle, n_cells_radial);
235 
236  // Create DG object, using max_poly = p+1 to allow for adjoint computation
237  std::shared_ptr < DGBase<dim, double> > dg = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree, poly_degree+1, grid);
238 
239  dg->allocate_system ();
240  // Initialize coarse grid solution with free-stream
241  dealii::VectorTools::interpolate(dg->dof_handler, initial_conditions, dg->solution);
242 
243  // Create ODE solver and ramp up the solution from p0
244  std::shared_ptr<ODE::ODESolverBase<dim, double>> ode_solver = ODE::ODESolverFactory<dim, double>::create_ODESolver(dg);
245  ode_solver->initialize_steady_polynomial_ramping(poly_degree);
246 
247  // setting up the target functional (error reduction)
248  std::shared_ptr< L2normError<dim, nstate, double> > L2normFunctional =
249  std::make_shared< L2normError<dim, nstate, double> >(dg,true,false);
250 
251  // initializing an adjoint for this case
252  Adjoint<dim, nstate, double> adjoint(dg, L2normFunctional, euler_physics_adtype);
253 
254  dealii::Vector<float> estimated_error_per_cell(grid->n_active_cells());
255  for (unsigned int igrid=0; igrid<n_grids; ++igrid) {
256 
257  // Interpolate solution from previous grid
258  if (igrid>0) {
259  dealii::LinearAlgebra::distributed::Vector<double> old_solution(dg->solution);
260  dealii::parallel::distributed::SolutionTransfer<dim, dealii::LinearAlgebra::distributed::Vector<double>, dealii::DoFHandler<dim>> solution_transfer(dg->dof_handler);
261  solution_transfer.prepare_for_coarsening_and_refinement(old_solution);
262  dg->high_order_grid->prepare_for_coarsening_and_refinement();
263 
264  // grid->refine_global (1);
265  // dealii::GridRefinement::refine_and_coarsen_fixed_fraction(*grid,
266  // dealii::parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(*grid,
267  dealii::parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(*grid,
268  // dealii::GridRefinement::refine_and_coarsen_fixed_number(*grid,
269  estimated_error_per_cell,
270  0.3,
271  0.03);
272  grid->execute_coarsening_and_refinement();
273  dg->high_order_grid->execute_coarsening_and_refinement();
274  dg->allocate_system ();
275  dg->solution.zero_out_ghosts();
276  solution_transfer.interpolate(dg->solution);
277  dg->solution.update_ghost_values();
278  }
279 
280  // std::string filename = "grid_cylinder-" + dealii::Utilities::int_to_string(igrid, 1) + ".eps";
281  // std::ofstream out (filename);
282  // dealii::GridOut grid_out;
283  // grid_out.write_eps (*grid, out);
284  // pcout << " Grid #" << igrid+1 << " . Number of cells: " << grid->n_global_active_cells() << std::endl;
285  // pcout << " written to " << filename << std::endl << std::endl;
286 
287 
288  const unsigned int n_global_active_cells = grid->n_global_active_cells();
289  const unsigned int n_dofs = dg->dof_handler.n_dofs();
290  pcout << "Dimension: " << dim << "\t Polynomial degree p: " << poly_degree << std::endl
291  << "Grid number: " << igrid+1 << "/" << n_grids
292  << ". Number of active cells: " << n_global_active_cells
293  << ". Number of degrees of freedom: " << n_dofs
294  << std::endl;
295 
296  //ode_solver->initialize_steady_polynomial_ramping (poly_degree);
297  // Solve the steady state problem
298  ode_solver->steady_state();
299 
300  const auto mapping = (*(dg->high_order_grid->mapping_fe_field));
301  dealii::hp::MappingCollection<dim> mapping_collection(mapping);
302  dealii::hp::FEValues<dim,dim> fe_values_collection_volume (mapping_collection, dg->fe_collection, dg->volume_quadrature_collection, dealii::update_values | dealii::update_JxW_values);
303  // Overintegrate the error to make sure there is not integration error in the error estimate
304  //int overintegrate = 0;
305  //dealii::QGauss<dim> quad_extra(dg->max_degree+1+overintegrate);
306  //dealii::FEValues<dim,dim> fe_values_extra(*(dg->high_order_grid->mapping_fe_field), dg->fe_collection[poly_degree], quad_extra,
307  // dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points);
308  std::array<double,nstate> soln_at_q;
309 
310  double l2error = 0;
311  double area = 0;
312  const double exact_area = (std::pow(outer_radius_adjoint+inner_radius_adjoint, 2.0) - std::pow(inner_radius_adjoint,2.0))*dealii::numbers::PI / 2.0;
313 
314 
315  const double entropy_inf = euler_physics_double->entropy_inf;
316 
317  estimated_error_per_cell.reinit(grid->n_active_cells());
318  // Integrate solution error and output error
319  for (auto cell = dg->dof_handler.begin_active(); cell!=dg->dof_handler.end(); ++cell) {
320 
321  if (!cell->is_locally_owned()) continue;
322 
323  const int i_fele = cell->active_fe_index();
324  const int i_quad = i_fele;
325  const int i_mapp = 0;
326 
327  fe_values_collection_volume.reinit (cell, i_quad, i_mapp, i_fele);
328  const dealii::FEValues<dim,dim> &fe_values_volume = fe_values_collection_volume.get_present_fe_values();
329 
330  //fe_values_volume.reinit (cell);
331  std::vector<dealii::types::global_dof_index> dofs_indices (fe_values_volume.dofs_per_cell);
332  cell->get_dof_indices (dofs_indices);
333 
334  double cell_l2error = 0;
335  for (unsigned int iquad=0; iquad<fe_values_volume.n_quadrature_points; ++iquad) {
336 
337  std::fill(soln_at_q.begin(), soln_at_q.end(), 0);
338  for (unsigned int idof=0; idof<fe_values_volume.dofs_per_cell; ++idof) {
339  const unsigned int istate = fe_values_volume.get_fe().system_to_component_index(idof).first;
340  soln_at_q[istate] += dg->solution[dofs_indices[idof]] * fe_values_volume.shape_value_component(idof, iquad, istate);
341  }
342  const double entropy = euler_physics_double->compute_entropy_measure(soln_at_q);
343 
344  const double uexact = entropy_inf;
345  cell_l2error += pow(entropy - uexact, 2) * fe_values_volume.JxW(iquad);
346 
347  area += fe_values_volume.JxW(iquad);
348  }
349 
350  // estimated_error_per_cell[cell->active_cell_index()] = cell_l2error;
351  l2error += cell_l2error;
352  }
353  const double l2error_mpi_sum = std::sqrt(dealii::Utilities::MPI::sum(l2error, mpi_communicator));
354  const double area_mpi_sum = dealii::Utilities::MPI::sum(area, mpi_communicator);
355 
356  // computing using Functional for comparison
357  const double l2error_functional = L2normFunctional->evaluate_functional(false,false);
358  pcout << "Error computed by original loop: " << l2error_mpi_sum << std::endl << "Error computed by the functional: " << std::sqrt(l2error_functional) << std::endl;
359 
360  // reinitializing the adjoint with the current values (from references)
361  adjoint.reinit();
362 
363  // evaluating the derivatives and the adjoint on the fine grid
364  adjoint.convert_to_state(PHiLiP::Adjoint<dim,nstate,double>::AdjointStateEnum::fine); // will do this automatically, but I prefer to repeat explicitly
365  adjoint.fine_grid_adjoint();
366  estimated_error_per_cell = adjoint.dual_weighted_residual(); // performing the error indicator computation
367 
368  // and outputing the fine properties
369  adjoint.output_results_vtk(igrid);
370 
371  adjoint.convert_to_state(PHiLiP::Adjoint<dim,nstate,double>::AdjointStateEnum::coarse); // this one is necessary though
372  adjoint.coarse_grid_adjoint();
373  adjoint.output_results_vtk(igrid);
374 
375  // Convergence table
376  const double dx = 1.0/pow(n_dofs,(1.0/dim));
377  grid_size[igrid] = dx;
378  entropy_error[igrid] = l2error_mpi_sum;
379 
380  convergence_table.add_value("p", poly_degree);
381  convergence_table.add_value("cells", n_global_active_cells);
382  convergence_table.add_value("DoFs", n_dofs);
383  convergence_table.add_value("dx", dx);
384  convergence_table.add_value("L2_entropy_error", l2error_mpi_sum);
385  convergence_table.add_value("area_error", std::abs(area_mpi_sum-exact_area));
386 
387 
388  pcout << " Grid size h: " << dx
389  << " L2-entropy_error: " << l2error_mpi_sum
390  << " Residual: " << ode_solver->residual_norm
391  << std::endl;
392 
393  if (igrid > 0) {
394  const double slope_soln_err = log(entropy_error[igrid]/entropy_error[igrid-1])
395  / log(grid_size[igrid]/grid_size[igrid-1]);
396  pcout << "From grid " << igrid-1
397  << " to grid " << igrid
398  << " dimension: " << dim
399  << " polynomial degree p: " << poly_degree
400  << std::endl
401  << " entropy_error1 " << entropy_error[igrid-1]
402  << " entropy_error2 " << entropy_error[igrid]
403  << " slope " << slope_soln_err
404  << std::endl;
405  }
406 
407  //output_results (igrid);
408  dg->output_results_vtk(igrid+10);
409  }
410 
411  pcout << " ********************************************" << std::endl
412  << " Convergence rates for p = " << poly_degree << std::endl
413  << " ********************************************" << std::endl;
414  convergence_table.evaluate_convergence_rates("L2_entropy_error", "cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
415  convergence_table.evaluate_convergence_rates("area_error", "cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
416  convergence_table.set_scientific("dx", true);
417  convergence_table.set_scientific("L2_entropy_error", true);
418  convergence_table.set_scientific("area_error", true);
419  if (pcout.is_active()) convergence_table.write_text(pcout.get_stream());
420 
421  convergence_table_vector.push_back(convergence_table);
422 
423  const double expected_slope = poly_degree+1;
424 
425  const double last_slope = log(entropy_error[n_grids-1]/entropy_error[n_grids-2])
426  / log(grid_size[n_grids-1]/grid_size[n_grids-2]);
427  double before_last_slope = last_slope;
428  if ( n_grids > 2 ) {
429  before_last_slope = log(entropy_error[n_grids-2]/entropy_error[n_grids-3])
430  / log(grid_size[n_grids-2]/grid_size[n_grids-3]);
431  }
432  const double slope_avg = 0.5*(before_last_slope+last_slope);
433  const double slope_diff = slope_avg-expected_slope;
434 
435  double slope_deficit_tolerance = -std::abs(manu_grid_conv_param.slope_deficit_tolerance);
436  if(poly_degree == 0) slope_deficit_tolerance *= 2; // Otherwise, grid sizes need to be much bigger for p=0
437 
438  if (slope_diff < slope_deficit_tolerance) {
439  pcout << std::endl
440  << "Convergence order not achieved. Average last 2 slopes of "
441  << slope_avg << " instead of expected "
442  << expected_slope << " within a tolerance of "
443  << slope_deficit_tolerance
444  << std::endl;
445  // p=0 just requires too many meshes to get into the asymptotic region.
446  if(poly_degree!=0) fail_conv_poly.push_back(poly_degree);
447  if(poly_degree!=0) fail_conv_slop.push_back(slope_avg);
448  }
449 
450  }
451  pcout << std::endl << std::endl << std::endl << std::endl;
452  pcout << " ********************************************" << std::endl;
453  pcout << " Convergence summary" << std::endl;
454  pcout << " ********************************************" << std::endl;
455  for (auto conv = convergence_table_vector.begin(); conv!=convergence_table_vector.end(); conv++) {
456  if (pcout.is_active()) conv->write_text(pcout.get_stream());
457  pcout << " ********************************************" << std::endl;
458  }
459  int n_fail_poly = fail_conv_poly.size();
460  if (n_fail_poly > 0) {
461  for (int ifail=0; ifail < n_fail_poly; ++ifail) {
462  const double expected_slope = fail_conv_poly[ifail]+1;
463  const double slope_deficit_tolerance = -0.1;
464  pcout << std::endl
465  << "Convergence order not achieved for polynomial p = "
466  << fail_conv_poly[ifail]
467  << ". Slope of "
468  << fail_conv_slop[ifail] << " instead of expected "
469  << expected_slope << " within a tolerance of "
470  << slope_deficit_tolerance
471  << std::endl;
472  }
473  }
474  return n_fail_poly;
475 }
476 
477 #if PHILIP_DIM==2
479 #endif
480 
481 } // Tests namespace
482 } // PHiLiP namespace
483 
484 
unsigned int degree_start
First polynomial degree to start the loop. If diffusion, must be at least 1.
PartialDifferentialEquation pde_type
Store the PDE type to be solved.
const bool uses_solution_gradient
Will evaluate solution gradient at quadrature points.
Definition: functional.h:315
real evaluate_volume_integrand(const PHiLiP::Physics::PhysicsBase< dim, nstate, real > &physics, const dealii::Point< dim, real > &phys_coord, const std::array< real, nstate > &soln_at_q, const std::array< dealii::Tensor< 1, dim, real >, nstate > &soln_grad_at_q) const override
non-template functions to override the template classes
const bool uses_solution_values
Will evaluate solution values at quadrature points.
Definition: functional.h:314
Function used to evaluate farfield conservative solution.
Parameters related to the manufactured convergence study.
Base class from which Advection, Diffusion, ConvectionDiffusion, and Euler is derived.
Definition: physics.h:34
Adjoint class.
Definition: adjoint.h:39
const MPI_Comm mpi_communicator
MPI communicator.
Definition: tests.h:39
L2normError(std::shared_ptr< PHiLiP::DGBase< dim, real >> dg_input, const bool uses_solution_values=true, const bool uses_solution_gradient=true)
Constructor.
Files for the baseline physics.
Definition: ADTypes.hpp:10
Sacado::Fad::DFad< real > FadType
Sacado AD type for first derivatives.
FadFadType evaluate_volume_integrand(const PHiLiP::Physics::PhysicsBase< dim, nstate, FadFadType > &physics, const dealii::Point< dim, FadFadType > &phys_coord, const std::array< FadFadType, nstate > &soln_at_q, const std::array< dealii::Tensor< 1, dim, FadFadType >, nstate > &soln_grad_at_q) const override
non-template functions to override the template classes
Sacado::Fad::DFad< FadType > FadFadType
Sacado AD type that allows 2nd derivatives.
int run_test() const override
Grid convergence on Euler Gaussian Bump.
real2 evaluate_volume_integrand(const PHiLiP::Physics::PhysicsBase< dim, nstate, real2 > &physics, const dealii::Point< dim, real2 > &, const std::array< real2, nstate > &soln_at_q, const std::array< dealii::Tensor< 1, dim, real2 >, nstate > &) const
Templated volume integrand of the functional, which is the point entropy generated squared...
unsigned int dimension
Number of dimensions. Note that it has to match the executable PHiLiP_xD.
Main parameter class that contains the various other sub-parameter classes.
ManufacturedConvergenceStudyParam manufactured_convergence_study_param
Contains parameters for manufactured convergence study.
const Parameters::AllParameters *const all_parameters
Pointer to all parameters.
Definition: tests.h:20
static std::shared_ptr< DGBase< dim, real, MeshType > > create_discontinuous_galerkin(const Parameters::AllParameters *const parameters_input, const unsigned int degree, const unsigned int max_degree_input, const unsigned int grid_degree_input, const std::shared_ptr< Triangulation > triangulation_input)
Creates a derived object DG, but returns it as DGBase.
Definition: dg_factory.cpp:10
const double entropy_inf
Entropy measure at infinity.
Definition: euler.h:127
Performs grid convergence for various polynomial degrees.
Euler equations. Derived from PhysicsBase.
Definition: euler.h:78
Create specified physics as PhysicsBase object.
static std::shared_ptr< PhysicsBase< dim, nstate, real > > create_Physics(const Parameters::AllParameters *const parameters_input, std::shared_ptr< ModelBase< dim, nstate, real > > model_input=nullptr)
Factory to return the correct physics given input file.
EulerCylinderAdjoint(const Parameters::AllParameters *const parameters_input)
Constructor.
dealii::ConditionalOStream pcout
ConditionalOStream.
Definition: tests.h:45
Functional base class.
Definition: functional.h:43
DGBase is independent of the number of state variables.
Definition: dg_base.hpp:82
std::vector< int > get_number_1d_cells(const int ngrids) const
Evaluates the number of cells to generate the grids for 1D grid based on input file.
Definition: tests.cpp:71
static std::shared_ptr< ODESolverBase< dim, real, MeshType > > create_ODESolver(std::shared_ptr< DGBase< dim, real, MeshType > > dg_input)
Creates either implicit or explicit ODE solver based on parameter value(no POD basis given) ...
Base class of all the tests.
Definition: tests.h:17
real compute_entropy_measure(const std::array< real, nstate > &conservative_soln) const
Evaluate entropy from conservative variables.
Definition: euler.cpp:288