5 #include <deal.II/base/convergence_table.h> 7 #include <deal.II/dofs/dof_tools.h> 9 #include <deal.II/grid/grid_generator.h> 10 #include <deal.II/grid/grid_refinement.h> 11 #include <deal.II/grid/grid_tools.h> 12 #include <deal.II/grid/grid_out.h> 13 #include <deal.II/grid/grid_in.h> 15 #include <deal.II/numerics/vector_tools.h> 16 #include <deal.II/numerics/data_out.h> 18 #include <deal.II/fe/fe_values.h> 19 #include <deal.II/fe/mapping_q.h> 21 #include "diffusion_exact_adjoint.h" 23 #include "parameters/all_parameters.h" 25 #include "physics/euler.h" 26 #include "physics/manufactured_solution.h" 27 #include "dg/dg_factory.hpp" 28 #include "dg/dg_base_state.hpp" 29 #include "ode_solver/ode_solver_factory.h" 31 #include "functional/functional.h" 32 #include "functional/adjoint.h" 34 #include "post_processor/physics_post_processor.h" 36 #include "ADTypes.hpp" 43 template <
int dim,
typename real>
50 val = (-1.0*pow(x,6)+3.0*pow(x,5)-3.0*pow(x,4)+pow(x,3));
52 real x = pos[0], y = pos[1];
53 val = (-1.0*pow(x,6)+3.0*pow(x,5)-3.0*pow(x,4)+pow(x,3))
54 * (-1.0*pow(y,6)+3.0*pow(y,5)-3.0*pow(y,4)+pow(y,3));
56 real x = pos[0], y = pos[1], z = pos[2];
57 val = (-1.0*pow(x,6)+3.0*pow(x,5)-3.0*pow(x,4)+pow(x,3))
58 * (-1.0*pow(y,6)+3.0*pow(y,5)-3.0*pow(y,4)+pow(y,3))
59 * (-1.0*pow(z,6)+3.0*pow(z,5)-3.0*pow(z,4)+pow(z,3));
66 template <
int dim,
typename real>
69 dealii::Tensor<1,dim,real> gradient;
73 gradient[0] = (-6.0*pow(x,5)+15.0*pow(x,4)-12.0*pow(x,3)+3.0*pow(x,2));
75 real x = pos[0], y = pos[1];
76 gradient[0] = (-6.0*pow(x,5)+15.0*pow(x,4)-12.0*pow(x,3)+3.0*pow(x,2))
77 * (-1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3));
78 gradient[1] = (-1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
79 * (-6.0*pow(y,5)+15.0*pow(y,4)-12.0*pow(y,3)+3.0*pow(y,2));
81 real x = pos[0], y = pos[1], z = pos[2];
82 gradient[0] = (-6.0*pow(x,5)+15.0*pow(x,4)-12.0*pow(x,3)+3.0*pow(x,2))
83 * (-1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
84 * (-1.0*pow(z,6)+ 3.0*pow(z,5)- 3.0*pow(z,4)+ pow(z,3));
85 gradient[1] = (-1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
86 * (-6.0*pow(y,5)+15.0*pow(y,4)-12.0*pow(y,3)+3.0*pow(y,2))
87 * (-1.0*pow(z,6)+ 3.0*pow(z,5)- 3.0*pow(z,4)+ pow(z,3));
88 gradient[2] = (-1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
89 * (-1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
90 * (-6.0*pow(z,5)+15.0*pow(z,4)-12.0*pow(z,3)+3.0*pow(z,2));
97 template <
int dim,
typename real>
100 const double pi = std::acos(-1);
108 real x = pos[0], y = pos[1];
112 real x = pos[0], y = pos[1], z = pos[2];
122 template <
int dim,
typename real>
125 const double pi = std::acos(-1);
127 dealii::Tensor<1,dim,real> gradient;
131 gradient[0] = pi*cos(pi*x);
133 real x = pos[0], y = pos[1];
134 gradient[0] = (pi*cos(pi*x))
136 gradient[1] = ( sin(pi*x))
139 real x = pos[0], y = pos[1], z = pos[2];
140 gradient[0] = (pi*cos(pi*x))
143 gradient[1] = ( sin(pi*x))
146 gradient[2] = ( sin(pi*x))
155 template <
int dim,
int nstate,
typename real>
157 const dealii::Point<dim,real> &pos,
158 const std::array<real,nstate> &,
160 const dealii::types::global_dof_index )
const 162 std::array<real,nstate> source;
168 val = (-30.0*pow(x,4)+60.0*pow(x,3)-36.0*pow(x,2)+6.0*x);
170 real x = pos[0], y = pos[1];
171 val = (-30.0*pow(x,4)+60.0*pow(x,3)-36.0*pow(x,2)+6.0*x)
172 * ( -1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
173 + ( -1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
174 * (-30.0*pow(y,4)+60.0*pow(y,3)-36.0*pow(y,2)+6.0*y);
176 real x = pos[0], y = pos[1], z = pos[2];
177 val = (-30.0*pow(x,4)+60.0*pow(x,3)-36.0*pow(x,2)+6.0*x)
178 * ( -1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
179 * ( -1.0*pow(z,6)+ 3.0*pow(z,5)- 3.0*pow(z,4)+ pow(z,3))
180 + ( -1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
181 * (-30.0*pow(y,4)+60.0*pow(y,3)-36.0*pow(y,2)+6.0*y)
182 * ( -1.0*pow(z,6)+ 3.0*pow(z,5)- 3.0*pow(z,4)+ pow(z,3))
183 + ( -1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
184 * ( -1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
185 * (-30.0*pow(z,4)+60.0*pow(z,3)-36.0*pow(z,2)+6.0*z);
188 for(
int istate = 0; istate < nstate; ++istate)
189 source[istate] = val;
194 template <
int dim,
int nstate,
typename real>
196 const dealii::Point<dim,real> &pos)
const 198 const double pi = std::acos(-1);
204 val = -pi*pi*sin(pi*x);
206 real x = pos[0], y = pos[1];
207 val = -pi*pi*sin(pi*x)
212 real x = pos[0], y = pos[1], z = pos[2];
213 val = -pi*pi*sin(pi*x)
227 template <
int dim,
int nstate,
typename real>
229 const dealii::Point<dim,real> &pos,
230 const std::array<real,nstate> &,
232 const dealii::types::global_dof_index )
const 234 const double pi = std::acos(-1);
236 std::array<real,nstate> source;
242 val = -pi*pi*sin(pi*x);
244 real x = pos[0], y = pos[1];
245 val = -pi*pi*sin(pi*x)
250 real x = pos[0], y = pos[1], z = pos[2];
251 val = -pi*pi*sin(pi*x)
262 for(
int istate = 0; istate < nstate; ++istate)
263 source[istate] = val;
268 template <
int dim,
int nstate,
typename real>
270 const dealii::Point<dim,real> &pos)
const 276 val = (-30.0*pow(x,4)+60.0*pow(x,3)-36.0*pow(x,2)+6.0*x);
278 real x = pos[0], y = pos[1];
279 val = (-30.0*pow(x,4)+60.0*pow(x,3)-36.0*pow(x,2)+6.0*x)
280 * ( -1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
281 + ( -1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
282 * (-30.0*pow(y,4)+60.0*pow(y,3)-36.0*pow(y,2)+6.0*y);
284 real x = pos[0], y = pos[1], z = pos[2];
285 val = (-30.0*pow(x,4)+60.0*pow(x,3)-36.0*pow(x,2)+6.0*x)
286 * ( -1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
287 * ( -1.0*pow(z,6)+ 3.0*pow(z,5)- 3.0*pow(z,4)+ pow(z,3))
288 + ( -1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
289 * (-30.0*pow(y,4)+60.0*pow(y,3)-36.0*pow(y,2)+6.0*y)
290 * ( -1.0*pow(z,6)+ 3.0*pow(z,5)- 3.0*pow(z,4)+ pow(z,3))
291 + ( -1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
292 * ( -1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
293 * (-30.0*pow(z,4)+60.0*pow(z,3)-36.0*pow(z,2)+6.0*z);
300 template <
int dim,
int nstate,
typename real>
301 template <
typename real2>
304 const dealii::Point<dim,real2> &phys_coord,
305 const std::array<real2,nstate> &soln_at_q,
306 const std::array<dealii::Tensor<1,dim,real2>,nstate> &)
const 317 for (
int istate=0; istate<nstate; ++istate) {
318 val += soln_at_q[istate] * objective_value;
325 template <
int dim,
int nstate>
329 template <
int dim,
int nstate>
332 pcout <<
"Running diffusion exact adjoint test case." << std::endl;
336 using GridEnum = ManParam::GridEnum;
346 const unsigned int p_start = manu_grid_conv_param.
degree_start;
347 const unsigned int p_end = manu_grid_conv_param.degree_end;
348 const unsigned int n_grids = manu_grid_conv_param.number_of_grids;
352 bool convection, diffusion;
353 if(pde_type == PdeEnum::diffusion){
357 pcout <<
"Can't run diffusion_exact_adjoint test case with other PDE types." << std::endl;
361 std::shared_ptr< Physics::PhysicsBase<dim, nstate, double> > physics_u_double
362 = std::make_shared< diffusion_u<dim, nstate, double> >(¶m, convection, diffusion);
363 std::shared_ptr< Physics::PhysicsBase<dim, nstate, double> > physics_v_double
364 = std::make_shared< diffusion_v<dim, nstate, double> >(¶m, convection, diffusion);
367 std::shared_ptr< Physics::PhysicsBase<dim, nstate, FadType> > physics_u_fadtype
368 = std::make_shared< diffusion_u<dim, nstate, FadType> >(¶m, convection, diffusion);
369 std::shared_ptr< Physics::PhysicsBase<dim, nstate, FadType> > physics_v_fadtype
370 = std::make_shared< diffusion_v<dim, nstate, FadType> >(¶m, convection, diffusion);
372 std::shared_ptr< Physics::PhysicsBase<dim, nstate, RadType> > physics_u_radtype
373 = std::make_shared< diffusion_u<dim, nstate, RadType> >(¶m, convection, diffusion);
374 std::shared_ptr< Physics::PhysicsBase<dim, nstate, RadType> > physics_v_radtype
375 = std::make_shared< diffusion_v<dim, nstate, RadType> >(¶m, convection, diffusion);
377 std::shared_ptr< Physics::PhysicsBase<dim, nstate, FadFadType> > physics_u_fadfadtype
378 = std::make_shared< diffusion_u<dim, nstate, FadFadType> >(¶m, convection, diffusion);
379 std::shared_ptr< Physics::PhysicsBase<dim, nstate, FadFadType> > physics_v_fadfadtype
380 = std::make_shared< diffusion_v<dim, nstate, FadFadType> >(¶m, convection, diffusion);
382 std::shared_ptr< Physics::PhysicsBase<dim, nstate, RadFadType> > physics_u_radfadtype
383 = std::make_shared< diffusion_u<dim, nstate, RadFadType> >(¶m, convection, diffusion);
384 std::shared_ptr< Physics::PhysicsBase<dim, nstate, RadFadType> > physics_v_radfadtype
385 = std::make_shared< diffusion_v<dim, nstate, RadFadType> >(¶m, convection, diffusion);
388 const double pi = std::acos(-1);
389 double exact_val = 0;
392 exact_val = (144*(pow(pi,2)-10)/pow(pi,5));
394 exact_val = 2 * (-144*(pow(pi,2)-10)/pow(pi,7)) * (144*(pow(pi,2)-10)/pow(pi,5));
396 exact_val = 3 * pow(-144*(pow(pi,2)-10)/pow(pi,7), 2) * (144*(pow(pi,2)-10)/pow(pi,5));
400 std::vector<dealii::ConvergenceTable> convergence_table_vector;
402 unsigned int n_fail = 0;
404 for(
unsigned int poly_degree = p_start; poly_degree <= p_end; ++poly_degree){
405 pcout <<
"Starting polynomial order: " << poly_degree << std::endl;
408 std::vector<double> grid_size(n_grids);
409 std::vector<double> output_error_u(n_grids);
410 std::vector<double> output_error_v(n_grids);
411 std::vector<double> soln_error_u(n_grids);
412 std::vector<double> soln_error_v(n_grids);
413 std::vector<double> adj_error_u(n_grids);
414 std::vector<double> adj_error_v(n_grids);
417 dealii::Vector<double> cellError_soln_u;
418 dealii::Vector<double> cellError_soln_v;
419 dealii::Vector<double> cellError_adj_u;
420 dealii::Vector<double> cellError_adj_v;
424 dealii::ConvergenceTable convergence_table;
427 using Triangulation = dealii::Triangulation<PHILIP_DIM>;
429 using Triangulation = dealii::parallel::distributed::Triangulation<PHILIP_DIM>;
431 std::shared_ptr <Triangulation> grid = std::make_shared<Triangulation> (
435 typename dealii::Triangulation<dim>::MeshSmoothing(
436 dealii::Triangulation<dim>::smoothing_on_refinement |
437 dealii::Triangulation<dim>::smoothing_on_coarsening));
440 const double left = 0.0;
441 const double right = 1.0;
443 for(
unsigned int igrid = 0; igrid < n_grids; ++igrid){
446 dealii::GridGenerator::subdivided_hyper_cube(*grid, n_1d_cells[igrid], left, right);
447 for (
auto cell = grid->begin_active(); cell != grid->end(); ++cell) {
448 cell->set_material_id(9002);
449 for (
unsigned int face=0; face<dealii::GeometryInfo<dim>::faces_per_cell; ++face) {
450 if (cell->face(face)->at_boundary()) cell->face(face)->set_boundary_id (1000);
464 dg_state_u->
set_physics(physics_u_double, physics_u_fadtype, physics_u_radtype, physics_u_fadfadtype, physics_u_radfadtype);
465 dg_state_v->set_physics(physics_v_double, physics_v_fadtype, physics_v_radtype, physics_v_fadfadtype, physics_v_radfadtype);
467 dg_u->allocate_system();
468 dg_v->allocate_system();
470 dg_u->solution *= 0.0;
471 dg_v->solution *= 0.0;
480 ode_solver_u->steady_state();
481 ode_solver_v->steady_state();
483 pcout <<
"Creating DiffusionFunctional... " << std::endl;
485 auto diffusion_functional_u = std::make_shared<DiffusionFunctional<dim,nstate,double>>(dg_u,physics_u_fadfadtype,
true,
false);
486 auto diffusion_functional_v = std::make_shared<DiffusionFunctional<dim,nstate,double>>(dg_v,physics_v_fadfadtype,
true,
false);
488 pcout <<
"Evaluating functional... " << std::endl;
490 double functional_val_u = diffusion_functional_u->evaluate_functional(
false,
false);
491 double functional_val_v = diffusion_functional_v->evaluate_functional(
false,
false);
494 pcout << std::endl <<
"Val1 = " << functional_val_u <<
"\tVal2 = " << functional_val_v << std::endl << std::endl;
497 double error_functional_u = std::abs(functional_val_u-exact_val);
498 double error_functional_v = std::abs(functional_val_v-exact_val);
500 pcout << std::endl <<
"error_val1 = " << error_functional_u <<
"\terror_val2 = " << error_functional_v << std::endl << std::endl;
507 pcout <<
"Solving for the discrete adjoints." << std::endl;
508 dealii::LinearAlgebra::distributed::Vector<double> adjoint_u = adj_u.
coarse_grid_adjoint();
509 dealii::LinearAlgebra::distributed::Vector<double> adjoint_v = adj_v.
coarse_grid_adjoint();
512 int overintegrate = 10;
513 dealii::QGauss<dim> quad_extra(poly_degree+overintegrate);
514 dealii::FEValues<dim,dim> fe_values_extra(*(dg_u->high_order_grid->mapping_fe_field), dg_u->fe_collection[poly_degree], quad_extra,
515 dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points);
516 const unsigned int n_quad_pts = fe_values_extra.n_quadrature_points;
519 double l2error_soln_u = 0;
520 double l2error_soln_v = 0;
523 double l2error_adj_u = 0;
524 double l2error_adj_v = 0;
527 std::array<double,nstate> soln_at_q_u;
528 std::array<double,nstate> soln_at_q_v;
529 std::array<double,nstate> adj_at_q_u;
530 std::array<double,nstate> adj_at_q_v;
533 cellError_soln_u.reinit(grid->n_active_cells());
534 cellError_soln_v.reinit(grid->n_active_cells());
535 cellError_adj_u.reinit(grid->n_active_cells());
536 cellError_adj_v.reinit(grid->n_active_cells());
538 std::vector<dealii::types::global_dof_index> dofs_indices(fe_values_extra.dofs_per_cell);
539 for(
auto cell = dg_u->dof_handler.begin_active(); cell != dg_u->dof_handler.end(); ++cell){
540 if(!cell->is_locally_owned())
continue;
542 fe_values_extra.reinit (cell);
543 cell->get_dof_indices(dofs_indices);
545 double cell_l2error_soln_u = 0;
546 double cell_l2error_soln_v = 0;
547 double cell_l2error_adj_u = 0;
548 double cell_l2error_adj_v = 0;
549 for(
unsigned int iquad = 0; iquad < n_quad_pts; ++iquad){
550 std::fill(soln_at_q_u.begin(), soln_at_q_u.end(), 0);
551 std::fill(soln_at_q_v.begin(), soln_at_q_v.end(), 0);
552 std::fill(adj_at_q_u.begin(), adj_at_q_u.end(), 0);
553 std::fill(adj_at_q_v.begin(), adj_at_q_v.end(), 0);
555 for (
unsigned int idof = 0; idof < fe_values_extra.dofs_per_cell; ++idof) {
556 const unsigned int istate = fe_values_extra.get_fe().system_to_component_index(idof).first;
557 soln_at_q_u[istate] += dg_u->solution[dofs_indices[idof]] * fe_values_extra.shape_value_component(idof, iquad, istate);
558 soln_at_q_v[istate] += dg_v->solution[dofs_indices[idof]] * fe_values_extra.shape_value_component(idof, iquad, istate);
559 adj_at_q_u[istate] += adjoint_u[dofs_indices[idof]] * fe_values_extra.shape_value_component(idof, iquad, istate);
560 adj_at_q_v[istate] += adjoint_v[dofs_indices[idof]] * fe_values_extra.shape_value_component(idof, iquad, istate);
563 const dealii::Point<dim,double> qpoint = (fe_values_extra.quadrature_point(iquad));
565 for (
int istate = 0; istate < nstate; ++istate){
566 const double soln_exact_u = physics_u_double->manufactured_solution_function->value(qpoint, istate);
567 const double soln_exact_v = physics_v_double->manufactured_solution_function->value(qpoint, istate);
570 cell_l2error_soln_u += pow(soln_at_q_u[istate] - soln_exact_u, 2) * fe_values_extra.JxW(iquad);
571 cell_l2error_soln_v += pow(soln_at_q_v[istate] - soln_exact_v, 2) * fe_values_extra.JxW(iquad);
574 cell_l2error_adj_u += pow(adj_at_q_u[istate] - soln_exact_v, 2) * fe_values_extra.JxW(iquad);
575 cell_l2error_adj_v += pow(adj_at_q_v[istate] - soln_exact_u, 2) * fe_values_extra.JxW(iquad);
584 cellError_soln_u[cell->active_cell_index()] = cell_l2error_soln_u;
585 cellError_soln_v[cell->active_cell_index()] = cell_l2error_soln_v;
586 cellError_adj_u[cell->active_cell_index()] = cell_l2error_adj_u;
587 cellError_adj_v[cell->active_cell_index()] = cell_l2error_adj_v;
590 l2error_soln_u += cell_l2error_soln_u;
591 l2error_soln_v += cell_l2error_soln_v;
592 l2error_adj_u += cell_l2error_adj_u;
593 l2error_adj_v += cell_l2error_adj_v;
595 const double l2error_soln_u_mpi_sum = std::sqrt(dealii::Utilities::MPI::sum(l2error_soln_u, mpi_communicator));
596 const double l2error_soln_v_mpi_sum = std::sqrt(dealii::Utilities::MPI::sum(l2error_soln_v, mpi_communicator));
597 const double l2error_adj_u_mpi_sum = std::sqrt(dealii::Utilities::MPI::sum(l2error_adj_u, mpi_communicator));
598 const double l2error_adj_v_mpi_sum = std::sqrt(dealii::Utilities::MPI::sum(l2error_adj_v, mpi_communicator));
602 const std::string filename_u =
"sol-u-" + std::to_string(igrid) +
".gnuplot";
603 std::ofstream gnuplot_output_u(filename_u);
604 dealii::DataOut<dim, dealii::DoFHandler<dim>> data_out_u;
605 data_out_u.attach_dof_handler(dg_u->dof_handler);
606 data_out_u.add_data_vector(dg_u->solution,
"u");
607 data_out_u.build_patches();
608 data_out_u.write_gnuplot(gnuplot_output_u);
610 const std::string filename_v =
"sol-v-" + std::to_string(igrid) +
".gnuplot";
611 std::ofstream gnuplot_output_v(filename_v);
612 dealii::DataOut<dim, dealii::DoFHandler<dim>> data_out_v;
613 data_out_v.attach_dof_handler(dg_v->dof_handler);
614 data_out_v.add_data_vector(dg_v->solution,
"u");
615 data_out_v.build_patches();
616 data_out_v.write_gnuplot(gnuplot_output_v);
618 pcout <<
"Outputting the results" << std::endl;
631 dealii::DataOut<dim, dealii::DoFHandler<dim>> data_out;
632 data_out.attach_dof_handler(dg_u->dof_handler);
640 dealii::Vector<float> subdomain(dg_u->triangulation->n_active_cells());
641 for (
unsigned int i = 0; i < subdomain.size(); ++i) {
642 subdomain(i) = dg_u->triangulation->locally_owned_subdomain();
644 data_out.add_data_vector(subdomain,
"subdomain", dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
646 std::vector<unsigned int> active_fe_indices;
647 dg_u->dof_handler.get_active_fe_indices(active_fe_indices);
648 dealii::Vector<double> active_fe_indices_dealiivector(active_fe_indices.begin(), active_fe_indices.end());
649 dealii::Vector<double> cell_poly_degree = active_fe_indices_dealiivector;
651 data_out.add_data_vector(active_fe_indices_dealiivector,
"PolynomialDegree", dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
653 std::vector<std::string> solution_names_u;
654 std::vector<std::string> solution_names_v;
655 std::vector<std::string> residual_names_u;
656 std::vector<std::string> residual_names_v;
657 std::vector<std::string> dIdw_names_u;
658 std::vector<std::string> dIdw_names_v;
659 std::vector<std::string> adjoint_names_u;
660 std::vector<std::string> adjoint_names_v;
661 for(
int s=0;s<nstate;++s) {
662 std::string varname0_u =
"state" + dealii::Utilities::int_to_string(s,1) +
"_u";
663 std::string varname0_v =
"state" + dealii::Utilities::int_to_string(s,1) +
"_v";
664 std::string varname1_u =
"residual" + dealii::Utilities::int_to_string(s,1) +
"_u";
665 std::string varname1_v =
"residual" + dealii::Utilities::int_to_string(s,1) +
"_v";
666 std::string varname2_u =
"dIdw" + dealii::Utilities::int_to_string(s,1) +
"_u";
667 std::string varname2_v =
"dIdw" + dealii::Utilities::int_to_string(s,1) +
"_v";
668 std::string varname3_u =
"psi" + dealii::Utilities::int_to_string(s,1) +
"_u";
669 std::string varname3_v =
"psi" + dealii::Utilities::int_to_string(s,1) +
"_v";
671 solution_names_u.push_back(varname0_u);
672 solution_names_v.push_back(varname0_v);
673 residual_names_u.push_back(varname1_u);
674 residual_names_v.push_back(varname1_v);
675 dIdw_names_u.push_back(varname2_u);
676 dIdw_names_v.push_back(varname2_v);
677 adjoint_names_u.push_back(varname3_u);
678 adjoint_names_v.push_back(varname3_v);
681 data_out.add_data_vector(dg_u->solution, solution_names_u, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
682 data_out.add_data_vector(dg_v->solution, solution_names_v, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
683 data_out.add_data_vector(dg_u->right_hand_side, residual_names_u, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
684 data_out.add_data_vector(dg_v->right_hand_side, residual_names_v, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
685 data_out.add_data_vector(adj_u.
dIdw_coarse, dIdw_names_u, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
686 data_out.add_data_vector(adj_v.
dIdw_coarse, dIdw_names_v, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
687 data_out.add_data_vector(adj_u.
adjoint_coarse, adjoint_names_u, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
688 data_out.add_data_vector(adj_v.
adjoint_coarse, adjoint_names_v, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
690 data_out.add_data_vector(cellError_soln_u,
"soln_u_err", dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
691 data_out.add_data_vector(cellError_soln_v,
"soln_v_err", dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
692 data_out.add_data_vector(cellError_adj_u,
"adj_u_err", dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
693 data_out.add_data_vector(cellError_adj_v,
"adj_v_err", dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
695 const int iproc = dealii::Utilities::MPI::this_mpi_process(mpi_communicator);
696 data_out.build_patches();
699 std::string filename =
"diffusion_exact_adjoint-" ;
700 filename += dealii::Utilities::int_to_string(dim, 1) +
"D-";
701 filename += dealii::Utilities::int_to_string(poly_degree, 1) +
"P-";
702 filename += dealii::Utilities::int_to_string(igrid, 4) +
".";
703 filename += dealii::Utilities::int_to_string(iproc, 4);
705 std::ofstream output(filename);
706 data_out.write_vtu(output);
709 std::vector<std::string> filenames;
710 for (
unsigned int iproc = 0; iproc < dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); ++iproc) {
711 std::string fn =
"diffusion_exact_adjoint-";
712 fn += dealii::Utilities::int_to_string(dim, 1) +
"D-";
713 fn += dealii::Utilities::int_to_string(poly_degree, 1) +
"P-";
714 fn += dealii::Utilities::int_to_string(igrid, 4) +
".";
715 fn += dealii::Utilities::int_to_string(iproc, 4);
717 filenames.push_back(fn);
719 std::string master_fn =
"diffusion_exact_adjoint-";
720 master_fn += dealii::Utilities::int_to_string(dim, 1) +
"D-";
721 master_fn += dealii::Utilities::int_to_string(poly_degree, 1) +
"P-";
722 master_fn += dealii::Utilities::int_to_string(igrid, 4) +
".pvtu";
723 std::ofstream master_output(master_fn);
724 data_out.write_pvtu_record(master_output, filenames);
729 const int n_dofs = dg_u->dof_handler.n_dofs();
730 const double dx = 1.0/pow(n_dofs,(1.0/dim));
733 convergence_table.add_value(
"p", poly_degree);
734 convergence_table.add_value(
"cells", grid->n_global_active_cells());
735 convergence_table.add_value(
"DoFs", n_dofs);
736 convergence_table.add_value(
"dx", dx);
737 convergence_table.add_value(
"soln_u_val", functional_val_u);
738 convergence_table.add_value(
"soln_v_val", functional_val_v);
739 convergence_table.add_value(
"soln_u_err", error_functional_u);
740 convergence_table.add_value(
"soln_v_err", error_functional_v);
741 convergence_table.add_value(
"soln_u_L2_err", l2error_soln_u_mpi_sum);
742 convergence_table.add_value(
"soln_v_L2_err", l2error_soln_v_mpi_sum);
743 convergence_table.add_value(
"adj_u_L2_err", l2error_adj_u_mpi_sum);
744 convergence_table.add_value(
"adj_v_L2_err", l2error_adj_v_mpi_sum);
747 grid_size[igrid] = dx;
748 output_error_u[igrid] = error_functional_u;
749 output_error_v[igrid] = error_functional_v;
750 soln_error_u[igrid] = l2error_soln_u_mpi_sum;
751 soln_error_v[igrid] = l2error_soln_v_mpi_sum;
752 adj_error_u[igrid] = l2error_adj_u_mpi_sum;
753 adj_error_v[igrid] = l2error_adj_v_mpi_sum;
757 convergence_table.evaluate_convergence_rates(
"soln_u_err",
"cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
758 convergence_table.evaluate_convergence_rates(
"soln_v_err",
"cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
759 convergence_table.evaluate_convergence_rates(
"soln_u_L2_err",
"cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
760 convergence_table.evaluate_convergence_rates(
"soln_v_L2_err",
"cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
761 convergence_table.evaluate_convergence_rates(
"adj_u_L2_err",
"cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
762 convergence_table.evaluate_convergence_rates(
"adj_v_L2_err",
"cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
763 convergence_table.set_scientific(
"dx",
true);
764 convergence_table.set_scientific(
"soln_u_val",
true);
765 convergence_table.set_scientific(
"soln_v_val",
true);
766 convergence_table.set_scientific(
"soln_u_err",
true);
767 convergence_table.set_scientific(
"soln_v_err",
true);
768 convergence_table.set_scientific(
"soln_u_L2_err",
true);
769 convergence_table.set_scientific(
"soln_v_L2_err",
true);
770 convergence_table.set_scientific(
"adj_u_L2_err",
true);
771 convergence_table.set_scientific(
"adj_v_L2_err",
true);
774 convergence_table_vector.push_back(convergence_table);
779 const double expected_slope_l2error_adj = poly_degree + 1;
786 const double avg_slope_error_l2error_adj_u = eval_avg_slope( adj_error_u, grid_size, n_grids);
787 const double avg_slope_error_l2error_adj_v = eval_avg_slope( adj_error_v, grid_size, n_grids);
794 const double diff_slope_error_l2error_adj_u = avg_slope_error_l2error_adj_u - expected_slope_l2error_adj;
795 const double diff_slope_error_l2error_adj_v = avg_slope_error_l2error_adj_v - expected_slope_l2error_adj;
798 double slope_deficit_tolerance = -std::abs(manu_grid_conv_param.slope_deficit_tolerance);
821 if(diff_slope_error_l2error_adj_u < slope_deficit_tolerance){
822 pcout <<
"Convergence order not achieved for l2error_adj_u." << std::endl
823 <<
"Average order of " << avg_slope_error_l2error_adj_u <<
" instead of expected " << expected_slope_l2error_adj << std::endl;
826 if(diff_slope_error_l2error_adj_v < slope_deficit_tolerance){
827 pcout <<
"Convergence order not achieved for l2error_adj_v." << std::endl
828 <<
"Average order of " << avg_slope_error_l2error_adj_v <<
" instead of expected " << expected_slope_l2error_adj << std::endl;
833 pcout << std::endl << std::endl << std::endl << std::endl;
834 pcout <<
" ********************************************" << std::endl;
835 pcout <<
" Convergence summary" << std::endl;
836 pcout <<
" ********************************************" << std::endl;
837 for (
auto conv = convergence_table_vector.begin(); conv!=convergence_table_vector.end(); conv++) {
838 if (
pcout.is_active()) conv->write_text(
pcout.get_stream());
839 pcout <<
" ********************************************" << std::endl;
846 double eval_avg_slope(std::vector<double> error, std::vector<double> grid_size,
unsigned int n_grids){
847 const double last_slope_error = log(error[n_grids-1]/error[n_grids-2])/(log(grid_size[n_grids-1]/grid_size[n_grids-2]));
848 double prev_slope_error = last_slope_error;
850 prev_slope_error = log(error[n_grids-2]/error[n_grids-3])/(log(grid_size[n_grids-2]/grid_size[n_grids-3]));
852 return 0.5*(last_slope_error+prev_slope_error);
real2 evaluate_volume_integrand(const PHiLiP::Physics::PhysicsBase< dim, nstate, real2 > &physics, const dealii::Point< dim, real2 > &phys_coord, const std::array< real2, nstate > &soln_at_q, const std::array< dealii::Tensor< 1, dim, real2 >, nstate > &soln_grad_at_q) const
Templated volume integrand.
real objective_function(const dealii::Point< dim, real > &pos) const override
objective function = f
unsigned int degree_start
First polynomial degree to start the loop. If diffusion, must be at least 1.
dealii::LinearAlgebra::distributed::Vector< real > dIdw_coarse
functional derivative (on the coarse grid)
PartialDifferentialEquation pde_type
Store the PDE type to be solved.
dealii::LinearAlgebra::distributed::Vector< real > adjoint_coarse
coarse grid adjoint ( )
virtual real objective_function(const dealii::Point< dim, real > &pos) const =0
defnined directly as part of the physics to make passing to the functional simpler ...
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.
PartialDifferentialEquation
Possible Partial Differential Equations to solve.
std::array< real, nstate > source_term(const dealii::Point< dim, real > &pos, const std::array< real, nstate > &solution, const real current_time, const dealii::types::global_dof_index cell_index) const override
source term = f
Files for the baseline physics.
DiffusionExactAdjoint(const Parameters::AllParameters *const parameters_input)
Constructor to call the TestsBase constructor to set parameters = parameters_input.
ManufacturedSolutionParam manufactured_solution_param
Associated manufactured solution parameters.
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.
real objective_function(const dealii::Point< dim, real > &pos) const override
objective function = g
ManufacturedConvergenceStudyParam manufactured_convergence_study_param
Contains parameters for manufactured convergence study.
dealii::LinearAlgebra::distributed::Vector< real > coarse_grid_adjoint()
Computes the coarse grid adjoint.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &pos, const unsigned int istate=0) const override
Gradient of the manufactured solution.
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.
real value(const dealii::Point< dim, real > &pos, const unsigned int istate=0) const override
overriding the function for the value and gradient
std::array< real, nstate > source_term(const dealii::Point< dim, real > &pos, const std::array< real, nstate > &solution, const real current_time, const dealii::types::global_dof_index cell_index) const override
source term = g
bool use_manufactured_source_term
Uses non-zero source term based on the manufactured solution and the PDE.
Abstract class templated on the number of state variables.
dealii::ConditionalOStream pcout
ConditionalOStream.
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &pos, const unsigned int istate=0) const override
Gradient of the manufactured solution.
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.
real value(const dealii::Point< dim, real > &pos, const unsigned int istate=0) const override
overriding the function for the value and gradient
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) ...
void set_physics(std::shared_ptr< Physics::PhysicsBase< dim, nstate, real > > pde_physics_double_input, std::shared_ptr< Physics::PhysicsBase< dim, nstate, FadType > > pde_physics_fad_input, std::shared_ptr< Physics::PhysicsBase< dim, nstate, RadType > > pde_physics_rad_input, std::shared_ptr< Physics::PhysicsBase< dim, nstate, FadFadType > > pde_physics_fad_fad_input, std::shared_ptr< Physics::PhysicsBase< dim, nstate, RadFadType > > pde_physics_rad_fad_input)
Base class of all the tests.
parent class to add the objective function directly to physics as a virtual class ...