1 #include "euler_cylinder_adjoint.h" 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> 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" 36 template <
int dim,
int nstate,
typename real>
49 template <
typename real2>
52 const dealii::Point<dim,real2> &,
53 const std::array<real2,nstate> &soln_at_q,
54 const std::array<dealii::Tensor<1,dim,real2>,nstate> &)
const 57 real2 cell_l2error = 0;
64 cell_l2error = pow(entropy - uexact, 2);
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 76 return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q);
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 89 return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q);
93 dealii::Point<2> center_adjoint(0.0,0.0);
94 const double inner_radius_adjoint = 1, outer_radius_adjoint = inner_radius_adjoint*20;
96 dealii::Point<2> warp_cylinder_adjoint (
const dealii::Point<2> &p)
98 const double rectangle_height = 1.0;
100 const double angle = p[1]/rectangle_height * dealii::numbers::PI;
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));
107 dealii::Point<2> q = p;
108 q[0] = -radius*cos(angle);
109 q[1] = radius*sin(angle);
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)
121 dealii::Point<2> p1(-1,0.0), p2(-0.0,1.0);
123 const bool colorize =
true;
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);
130 dealii::GridTools::transform (&warp_cylinder_adjoint, tria);
132 tria.set_all_manifold_ids(0);
133 tria.set_manifold(0, dealii::SphericalManifold<2>(center_adjoint));
136 for (
auto cell = tria.begin_active(); cell != tria.end(); ++cell) {
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);
143 }
else if (current_id == 1) {
144 cell->face(face)->set_boundary_id (1001);
145 }
else if (current_id == 2) {
146 cell->face(face)->set_boundary_id (1001);
147 }
else if (current_id == 3) {
148 cell->face(face)->set_boundary_id (1001);
157 template <
int dim,
int nstate>
163 template<
int dim,
int nstate>
167 pcout <<
" Running Euler cylinder adjoint entropy convergence. " << std::endl;
168 static_assert(dim==2);
170 using GridEnum = ManParam::GridEnum;
174 Assert(param.
pde_type == param.PartialDifferentialEquation::euler, dealii::ExcMessage(
"Can't run Euler case if PDE is not Euler"));
179 const unsigned int p_start = manu_grid_conv_param.
degree_start;
180 const unsigned int p_end = manu_grid_conv_param.degree_end;
182 const unsigned int n_grids_input = manu_grid_conv_param.number_of_grids;
184 std::shared_ptr< Physics::PhysicsBase<dim,nstate,double> > physics_double
187 std::shared_ptr< Physics::Euler<dim,nstate,double> > euler_physics_double
190 std::shared_ptr< Physics::PhysicsBase<dim,nstate,Sacado::Fad::DFad<double>> > euler_physics_adtype
211 std::vector<int> fail_conv_poly;
212 std::vector<double> fail_conv_slop;
213 std::vector<dealii::ConvergenceTable> convergence_table_vector;
215 for (
unsigned int poly_degree = p_start; poly_degree <= p_end; ++poly_degree) {
218 unsigned int n_grids = n_grids_input;
219 if (poly_degree <= 1) n_grids = n_grids_input;
221 std::vector<double> entropy_error(n_grids);
222 std::vector<double> grid_size(n_grids);
226 dealii::ConvergenceTable convergence_table;
229 using Triangulation = dealii::parallel::distributed::Triangulation<dim>;
230 std::shared_ptr <Triangulation> grid = std::make_shared<Triangulation> (this->
mpi_communicator);
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);
239 dg->allocate_system ();
241 dealii::VectorTools::interpolate(dg->dof_handler, initial_conditions, dg->solution);
245 ode_solver->initialize_steady_polynomial_ramping(poly_degree);
248 std::shared_ptr< L2normError<dim, nstate, double> > L2normFunctional =
249 std::make_shared< L2normError<dim, nstate, double> >(dg,
true,
false);
254 dealii::Vector<float> estimated_error_per_cell(grid->n_active_cells());
255 for (
unsigned int igrid=0; igrid<n_grids; ++igrid) {
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();
267 dealii::parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(*grid,
269 estimated_error_per_cell,
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();
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
298 ode_solver->steady_state();
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);
308 std::array<double,nstate> soln_at_q;
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;
315 const double entropy_inf = euler_physics_double->entropy_inf;
317 estimated_error_per_cell.reinit(grid->n_active_cells());
319 for (
auto cell = dg->dof_handler.begin_active(); cell!=dg->dof_handler.end(); ++cell) {
321 if (!cell->is_locally_owned())
continue;
323 const int i_fele = cell->active_fe_index();
324 const int i_quad = i_fele;
325 const int i_mapp = 0;
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();
331 std::vector<dealii::types::global_dof_index> dofs_indices (fe_values_volume.dofs_per_cell);
332 cell->get_dof_indices (dofs_indices);
334 double cell_l2error = 0;
335 for (
unsigned int iquad=0; iquad<fe_values_volume.n_quadrature_points; ++iquad) {
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);
342 const double entropy = euler_physics_double->compute_entropy_measure(soln_at_q);
344 const double uexact = entropy_inf;
345 cell_l2error += pow(entropy - uexact, 2) * fe_values_volume.JxW(iquad);
347 area += fe_values_volume.JxW(iquad);
351 l2error += cell_l2error;
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);
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;
365 adjoint.fine_grid_adjoint();
366 estimated_error_per_cell = adjoint.dual_weighted_residual();
369 adjoint.output_results_vtk(igrid);
372 adjoint.coarse_grid_adjoint();
373 adjoint.output_results_vtk(igrid);
376 const double dx = 1.0/pow(n_dofs,(1.0/dim));
377 grid_size[igrid] = dx;
378 entropy_error[igrid] = l2error_mpi_sum;
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));
388 pcout <<
" Grid size h: " << dx
389 <<
" L2-entropy_error: " << l2error_mpi_sum
390 <<
" Residual: " << ode_solver->residual_norm
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
401 <<
" entropy_error1 " << entropy_error[igrid-1]
402 <<
" entropy_error2 " << entropy_error[igrid]
403 <<
" slope " << slope_soln_err
408 dg->output_results_vtk(igrid+10);
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());
421 convergence_table_vector.push_back(convergence_table);
423 const double expected_slope = poly_degree+1;
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;
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]);
432 const double slope_avg = 0.5*(before_last_slope+last_slope);
433 const double slope_diff = slope_avg-expected_slope;
435 double slope_deficit_tolerance = -std::abs(manu_grid_conv_param.slope_deficit_tolerance);
436 if(poly_degree == 0) slope_deficit_tolerance *= 2;
438 if (slope_diff < slope_deficit_tolerance) {
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
446 if(poly_degree!=0) fail_conv_poly.push_back(poly_degree);
447 if(poly_degree!=0) fail_conv_slop.push_back(slope_avg);
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;
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;
465 <<
"Convergence order not achieved for polynomial p = " 466 << fail_conv_poly[ifail]
468 << fail_conv_slop[ifail] <<
" instead of expected " 469 << expected_slope <<
" within a tolerance of " 470 << slope_deficit_tolerance
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.
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.
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.
const MPI_Comm mpi_communicator
MPI communicator.
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.
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.
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.
const double entropy_inf
Entropy measure at infinity.
Performs grid convergence for various polynomial degrees.
Euler equations. Derived from PhysicsBase.
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.
DGBase is independent of the number of state variables.
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.
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.
real compute_entropy_measure(const std::array< real, nstate > &conservative_soln) const
Evaluate entropy from conservative variables.