2 #include "functional.h" 4 #include <deal.II/base/function.h> 5 #include <deal.II/base/symmetric_tensor.h> 6 #include <deal.II/differentiation/ad/sacado_math.h> 7 #include <deal.II/differentiation/ad/sacado_number_types.h> 8 #include <deal.II/differentiation/ad/sacado_product_types.h> 9 #include <deal.II/dofs/dof_tools.h> 10 #include <deal.II/fe/fe_q.h> 11 #include <deal.II/fe/fe_values.h> 12 #include <deal.II/lac/la_parallel_vector.h> 18 #include "dg/dg_base_state.hpp" 19 #include "lift_drag.hpp" 20 #include "physics/model.h" 21 #include "physics/model_factory.h" 22 #include "physics/physics.h" 23 #include "physics/physics_factory.h" 30 template<
int dim,
typename real1,
typename real2>
31 dealii::Tensor<1,dim,real1> vmult(
const dealii::Tensor<2,dim,real1> A,
const dealii::Tensor<1,dim,real2> x)
33 dealii::Tensor<1,dim,real1> y;
34 for (
int row=0;row<dim;++row) {
36 for (
int col=0;col<dim;++col) {
37 y[row] += A[row][col] * x[col];
48 template<
int dim,
typename real1>
49 real1 norm(
const dealii::Tensor<1,dim,real1> x)
52 for (
int row=0;row<dim;++row) {
53 val += x[row] * x[row];
60 template <
int dim,
int nstate,
typename real,
typename MeshType>
61 FunctionalNormLpVolume<dim,nstate,real,MeshType>::FunctionalNormLpVolume(
63 std::shared_ptr<DGBase<dim,real,MeshType>> _dg,
64 const bool _uses_solution_values,
65 const bool _uses_solution_gradient) :
66 Functional<dim,nstate,real,MeshType>::Functional(_dg, _uses_solution_values, _uses_solution_gradient),
69 template <
int dim,
int nstate,
typename real,
typename MeshType>
70 template <
typename real2>
71 real2 FunctionalNormLpVolume<dim,nstate,real,MeshType>::evaluate_volume_integrand(
73 const dealii::Point<dim,real2> & ,
74 const std::array<real2,nstate> & soln_at_q,
75 const std::array<dealii::Tensor<1,dim,real2>,nstate> &)
const 77 real2 lpnorm_value = 0;
78 for(
unsigned int istate = 0; istate < nstate; ++istate)
79 lpnorm_value += pow(abs(soln_at_q[istate]), this->normLp);
83 template <
int dim,
int nstate,
typename real,
typename MeshType>
84 FunctionalNormLpBoundary<dim,nstate,real,MeshType>::FunctionalNormLpBoundary(
86 std::vector<unsigned int> _boundary_vector,
87 const bool _use_all_boundaries,
88 std::shared_ptr<DGBase<dim,real,MeshType>> _dg,
89 const bool _uses_solution_values,
90 const bool _uses_solution_gradient) :
91 Functional<dim,nstate,real,MeshType>::Functional(_dg, _uses_solution_values, _uses_solution_gradient),
93 boundary_vector(_boundary_vector),
94 use_all_boundaries(_use_all_boundaries) {}
96 template <
int dim,
int nstate,
typename real,
typename MeshType>
97 template <
typename real2>
98 real2 FunctionalNormLpBoundary<dim,nstate,real,MeshType>::evaluate_boundary_integrand(
100 const unsigned int boundary_id,
101 const dealii::Point<dim,real2> & ,
102 const dealii::Tensor<1,dim,real2> & ,
103 const std::array<real2,nstate> & soln_at_q,
104 const std::array<dealii::Tensor<1,dim,real2>,nstate> &)
const 106 real2 lpnorm_value = 0;
109 auto boundary_vector_index = std::find(this->boundary_vector.begin(), this->boundary_vector.end(), boundary_id);
110 bool eval_boundary = this->use_all_boundaries || boundary_vector_index != this->boundary_vector.end();
115 for(
unsigned int istate = 0; istate < nstate; ++istate)
116 lpnorm_value += pow(abs(soln_at_q[istate]), this->normLp);
121 template <
int dim,
int nstate,
typename real,
typename MeshType>
122 FunctionalWeightedIntegralVolume<dim,nstate,real,MeshType>::FunctionalWeightedIntegralVolume(
123 std::shared_ptr<ManufacturedSolutionFunction<dim,real>> _weight_function_double,
124 std::shared_ptr<ManufacturedSolutionFunction<dim,FadFadType>> _weight_function_adtype,
125 const bool _use_weight_function_laplacian,
126 std::shared_ptr<DGBase<dim,real,MeshType>> _dg,
127 const bool _uses_solution_values,
128 const bool _uses_solution_gradient) :
129 Functional<dim,nstate,real,MeshType>::Functional(_dg, _uses_solution_values, _uses_solution_gradient),
130 weight_function_double(_weight_function_double),
131 weight_function_adtype(_weight_function_adtype),
132 use_weight_function_laplacian(_use_weight_function_laplacian) {}
134 template <
int dim,
int nstate,
typename real,
typename MeshType>
135 template <
typename real2>
136 real2 FunctionalWeightedIntegralVolume<dim,nstate,real,MeshType>::evaluate_volume_integrand(
138 const dealii::Point<dim,real2> & phys_coord,
139 const std::array<real2,nstate> & soln_at_q,
140 const std::array<dealii::Tensor<1,dim,real2>,nstate> & ,
141 std::shared_ptr<ManufacturedSolutionFunction<dim,real2>> weight_function)
const 145 if(this->use_weight_function_laplacian){
146 for(
unsigned int istate = 0; istate < nstate; ++istate)
147 val += soln_at_q[istate] * dealii::trace(weight_function->hessian(phys_coord, istate));
149 for(
unsigned int istate = 0; istate < nstate; ++istate)
150 val += soln_at_q[istate] * weight_function->value(phys_coord, istate);
156 template <
int dim,
int nstate,
typename real,
typename MeshType>
157 FunctionalWeightedIntegralBoundary<dim,nstate,real,MeshType>::FunctionalWeightedIntegralBoundary(
158 std::shared_ptr<ManufacturedSolutionFunction<dim,real>> _weight_function_double,
159 std::shared_ptr<ManufacturedSolutionFunction<dim,FadFadType>> _weight_function_adtype,
160 const bool _use_weight_function_laplacian,
161 std::vector<unsigned int> _boundary_vector,
162 const bool _use_all_boundaries,
163 std::shared_ptr<DGBase<dim,real,MeshType>> _dg,
164 const bool _uses_solution_values,
165 const bool _uses_solution_gradient) :
166 Functional<dim,nstate,real,MeshType>::Functional(_dg, _uses_solution_values, _uses_solution_gradient),
167 weight_function_double(_weight_function_double),
168 weight_function_adtype(_weight_function_adtype),
169 use_weight_function_laplacian(_use_weight_function_laplacian),
170 boundary_vector(_boundary_vector),
171 use_all_boundaries(_use_all_boundaries) {}
173 template <
int dim,
int nstate,
typename real,
typename MeshType>
174 template <
typename real2>
175 real2 FunctionalWeightedIntegralBoundary<dim,nstate,real,MeshType>::evaluate_boundary_integrand(
177 const unsigned int boundary_id,
178 const dealii::Point<dim,real2> & phys_coord,
179 const dealii::Tensor<1,dim,real2> & ,
180 const std::array<real2,nstate> & soln_at_q,
181 const std::array<dealii::Tensor<1,dim,real2>,nstate> &,
182 std::shared_ptr<ManufacturedSolutionFunction<dim,real2>> weight_function)
const 187 auto boundary_vector_index = std::find(this->boundary_vector.begin(), this->boundary_vector.end(), boundary_id);
188 bool eval_boundary = this->use_all_boundaries || boundary_vector_index != this->boundary_vector.end();
193 if(this->use_weight_function_laplacian){
194 for(
unsigned int istate = 0; istate < nstate; ++istate)
195 val += soln_at_q[istate] * dealii::trace(weight_function->hessian(phys_coord, istate));
197 for(
unsigned int istate = 0; istate < nstate; ++istate)
198 val += soln_at_q[istate] * weight_function->value(phys_coord, istate);
204 template <
int dim,
int nstate,
typename real,
typename MeshType>
205 FunctionalErrorNormLpVolume<dim,nstate,real,MeshType>::FunctionalErrorNormLpVolume(
206 const double _normLp,
207 std::shared_ptr<DGBase<dim,real,MeshType>> _dg,
208 const bool _uses_solution_values,
209 const bool _uses_solution_gradient) :
210 Functional<dim,nstate,real,MeshType>::Functional(_dg, _uses_solution_values, _uses_solution_gradient),
213 template <
int dim,
int nstate,
typename real,
typename MeshType>
214 template <
typename real2>
215 real2 FunctionalErrorNormLpVolume<dim,nstate,real,MeshType>::evaluate_volume_integrand(
217 const dealii::Point<dim,real2> & phys_coord,
218 const std::array<real2,nstate> & soln_at_q,
219 const std::array<dealii::Tensor<1,dim,real2>,nstate> &)
const 221 real2 lpnorm_value = 0;
222 for(
unsigned int istate = 0; istate < nstate; ++istate){
224 lpnorm_value += pow(abs(soln_at_q[istate] - uexact), this->normLp);
229 template <
int dim,
int nstate,
typename real,
typename MeshType>
230 FunctionalErrorNormLpBoundary<dim,nstate,real,MeshType>::FunctionalErrorNormLpBoundary(
231 const double _normLp,
232 std::vector<unsigned int> _boundary_vector,
233 const bool _use_all_boundaries,
234 std::shared_ptr<DGBase<dim,real,MeshType>> _dg,
235 const bool _uses_solution_values,
236 const bool _uses_solution_gradient) :
237 Functional<dim,nstate,real,MeshType>::Functional(_dg, _uses_solution_values, _uses_solution_gradient),
239 boundary_vector(_boundary_vector),
240 use_all_boundaries(_use_all_boundaries) {}
242 template <
int dim,
int nstate,
typename real,
typename MeshType>
243 template <
typename real2>
244 real2 FunctionalErrorNormLpBoundary<dim,nstate,real,MeshType>::evaluate_boundary_integrand(
246 const unsigned int boundary_id,
247 const dealii::Point<dim,real2> & phys_coord,
248 const dealii::Tensor<1,dim,real2> & ,
249 const std::array<real2,nstate> & soln_at_q,
250 const std::array<dealii::Tensor<1,dim,real2>,nstate> &)
const 252 real2 lpnorm_value = 0;
255 auto boundary_vector_index = std::find(this->boundary_vector.begin(), this->boundary_vector.end(), boundary_id);
256 bool eval_boundary = this->use_all_boundaries || boundary_vector_index != this->boundary_vector.end();
261 for(
int istate = 0; istate < nstate; ++istate){
263 lpnorm_value += pow(abs(soln_at_q[istate] - uexact), this->normLp);
270 template <
int dim,
int nstate,
typename real,
typename MeshType>
271 template <
typename real2>
274 const dealii::Point<dim,real2> &,
275 const std::array<real2,nstate> &soln_at_q,
276 const std::array<dealii::Tensor<1,dim,real2>,nstate> &)
const 281 for (
int istate=0; istate<nstate; ++istate) {
282 val += soln_at_q[istate];
288 template <
int dim,
int nstate,
typename real,
typename MeshType>
291 const bool uses_solution_values,
292 const bool uses_solution_gradient)
296 template <
int dim,
int nstate,
typename real,
typename MeshType>
297 template <
typename real2>
300 const unsigned int boundary_id,
301 const dealii::Point<dim,real2> & ,
302 const dealii::Tensor<1,dim,real2> & ,
303 const std::array<real2,nstate> & soln_at_q,
304 const std::array<dealii::Tensor<1,dim,real2>,nstate> & )
const 308 if(boundary_id == 1002){
320 template <
int dim,
int nstate,
typename real,
typename MeshType>
323 const bool _uses_solution_values,
324 const bool _uses_solution_gradient)
326 , d2IdWdW(std::make_shared<dealii::TrilinosWrappers::SparseMatrix>())
327 , d2IdWdX(std::make_shared<dealii::TrilinosWrappers::SparseMatrix>())
328 , d2IdXdX(std::make_shared<dealii::TrilinosWrappers::SparseMatrix>())
329 , uses_solution_values(_uses_solution_values)
330 , uses_solution_gradient(_uses_solution_gradient)
331 ,
pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0)
333 using FadType = Sacado::Fad::DFad<real>;
334 using FadFadType = Sacado::Fad::DFad<FadType>;
340 template <
int dim,
int nstate,
typename real,
typename MeshType>
343 solution_value.reinit(
dg->solution);
344 solution_value *= 0.0;
345 volume_nodes_value.reinit(
dg->high_order_grid->volume_nodes);
346 volume_nodes_value *= 0.0;
348 solution_dIdW.reinit(
dg->solution);
349 solution_dIdW *= 0.0;
350 volume_nodes_dIdW.reinit(
dg->high_order_grid->volume_nodes);
351 volume_nodes_dIdW *= 0.0;
353 solution_dIdX.reinit(
dg->solution);
354 solution_dIdX *= 0.0;
355 volume_nodes_dIdX.reinit(
dg->high_order_grid->volume_nodes);
356 volume_nodes_dIdX *= 0.0;
358 solution_d2I.reinit(
dg->solution);
360 volume_nodes_d2I.reinit(
dg->high_order_grid->volume_nodes);
361 volume_nodes_d2I *= 0.0;
364 template <
int dim,
int nstate,
typename real,
typename MeshType>
368 const bool _uses_solution_values,
369 const bool _uses_solution_gradient)
370 :
Functional(_dg, _uses_solution_values, _uses_solution_gradient)
372 physics_fad_fad = _physics_fad_fad;
375 template <
int dim,
int nstate,
typename real,
typename MeshType>
378 dg->solution = solution_set;
381 template <
int dim,
int nstate,
typename real,
typename MeshType>
384 dg->high_order_grid->volume_nodes = volume_nodes_set;
387 template <
int dim,
int nstate,
typename real,
typename MeshType>
391 dealii::IndexSet locally_owned_dofs = dg->high_order_grid->dof_handler_grid.locally_owned_dofs();
392 dealii::IndexSet locally_relevant_dofs, ghost_dofs;
393 dealii::DoFTools::extract_locally_relevant_dofs(dg->high_order_grid->dof_handler_grid, locally_relevant_dofs);
394 ghost_dofs = locally_relevant_dofs;
395 ghost_dofs.subtract_set(locally_owned_dofs);
396 dIdX.reinit(locally_owned_dofs, ghost_dofs, MPI_COMM_WORLD);
399 template <
int dim,
int nstate,
typename real,
typename MeshType>
404 dealii::IndexSet locally_owned_dofs = dg->dof_handler.locally_owned_dofs();
405 dIdw.reinit(locally_owned_dofs, MPI_COMM_WORLD);
412 dealii::SparsityPattern sparsity_pattern_d2IdWdX = dg->get_d2RdWdX_sparsity_pattern ();
413 const dealii::IndexSet &row_parallel_partitioning_d2IdWdX = dg->locally_owned_dofs;
414 const dealii::IndexSet &col_parallel_partitioning_d2IdWdX = dg->high_order_grid->locally_owned_dofs_grid;
415 d2IdWdX->reinit(row_parallel_partitioning_d2IdWdX, col_parallel_partitioning_d2IdWdX, sparsity_pattern_d2IdWdX, MPI_COMM_WORLD);
419 dealii::SparsityPattern sparsity_pattern_d2IdWdW = dg->get_d2RdWdW_sparsity_pattern ();
420 const dealii::IndexSet &row_parallel_partitioning_d2IdWdW = dg->locally_owned_dofs;
421 const dealii::IndexSet &col_parallel_partitioning_d2IdWdW = dg->locally_owned_dofs;
422 d2IdWdW->reinit(row_parallel_partitioning_d2IdWdW, col_parallel_partitioning_d2IdWdW, sparsity_pattern_d2IdWdW, MPI_COMM_WORLD);
426 dealii::SparsityPattern sparsity_pattern_d2IdXdX = dg->get_d2RdXdX_sparsity_pattern ();
427 const dealii::IndexSet &row_parallel_partitioning_d2IdXdX = dg->high_order_grid->locally_owned_dofs_grid;
428 const dealii::IndexSet &col_parallel_partitioning_d2IdXdX = dg->high_order_grid->locally_owned_dofs_grid;
429 d2IdXdX->reinit(row_parallel_partitioning_d2IdXdX, col_parallel_partitioning_d2IdXdX, sparsity_pattern_d2IdXdX, MPI_COMM_WORLD);
435 template <
int dim,
int nstate,
typename real,
typename MeshType>
437 const bool compute_dIdW,
const bool compute_dIdX,
const bool compute_d2I,
438 const Sacado::Fad::DFad<Sacado::Fad::DFad<real>> volume_local_sum,
439 std::vector<dealii::types::global_dof_index> cell_soln_dofs_indices,
440 std::vector<dealii::types::global_dof_index> cell_metric_dofs_indices)
442 using FadType = Sacado::Fad::DFad<real>;
444 const unsigned int n_total_indep = volume_local_sum.size();
445 (void) n_total_indep;
446 const unsigned int n_soln_dofs_cell = cell_soln_dofs_indices.size();
447 const unsigned int n_metric_dofs_cell = cell_metric_dofs_indices.size();
448 unsigned int i_derivative = 0;
451 std::vector<real> local_dIdw(n_soln_dofs_cell);
452 for(
unsigned int idof = 0; idof < n_soln_dofs_cell; ++idof){
453 local_dIdw[idof] = volume_local_sum.dx(i_derivative++).val();
455 dIdw.add(cell_soln_dofs_indices, local_dIdw);
458 std::vector<real> local_dIdX(n_metric_dofs_cell);
459 for(
unsigned int idof = 0; idof < n_metric_dofs_cell; ++idof){
460 local_dIdX[idof] = volume_local_sum.dx(i_derivative++).val();
462 dIdX.add(cell_metric_dofs_indices, local_dIdX);
464 if (compute_dIdW || compute_dIdX) AssertDimension(i_derivative, n_total_indep);
466 std::vector<real> dWidW(n_soln_dofs_cell);
467 std::vector<real> dWidX(n_metric_dofs_cell);
468 std::vector<real> dXidX(n_metric_dofs_cell);
472 for (
unsigned int idof=0; idof<n_soln_dofs_cell; ++idof) {
474 unsigned int j_derivative = 0;
475 const FadType dWi = volume_local_sum.dx(i_derivative++);
477 for (
unsigned int jdof=0; jdof<n_soln_dofs_cell; ++jdof) {
478 dWidW[jdof] = dWi.dx(j_derivative++);
480 d2IdWdW->add(cell_soln_dofs_indices[idof], cell_soln_dofs_indices, dWidW);
482 for (
unsigned int jdof=0; jdof<n_metric_dofs_cell; ++jdof) {
483 dWidX[jdof] = dWi.dx(j_derivative++);
485 d2IdWdX->add(cell_soln_dofs_indices[idof], cell_metric_dofs_indices, dWidX);
488 for (
unsigned int idof=0; idof<n_metric_dofs_cell; ++idof) {
490 const FadType dXi = volume_local_sum.dx(i_derivative++);
492 unsigned int j_derivative = n_soln_dofs_cell;
493 for (
unsigned int jdof=0; jdof<n_metric_dofs_cell; ++jdof) {
494 dXidX[jdof] = dXi.dx(j_derivative++);
496 d2IdXdX->add(cell_metric_dofs_indices[idof], cell_metric_dofs_indices, dXidX);
499 AssertDimension(i_derivative, n_total_indep);
502 template <
int dim,
int nstate,
typename real,
typename MeshType>
503 template <
typename real2>
506 const std::vector< real2 > &soln_coeff,
507 const dealii::FESystem<dim> &fe_solution,
508 const std::vector< real2 > &coords_coeff,
509 const dealii::FESystem<dim> &fe_metric,
510 const dealii::Quadrature<dim> &volume_quadrature)
const 512 const unsigned int n_vol_quad_pts = volume_quadrature.size();
513 const unsigned int n_soln_dofs_cell = soln_coeff.size();
514 const unsigned int n_metric_dofs_cell = coords_coeff.size();
516 real2 volume_local_sum = 0.0;
517 for (
unsigned int iquad=0; iquad<n_vol_quad_pts; ++iquad) {
519 const dealii::Point<dim,double> &ref_point = volume_quadrature.point(iquad);
520 const double quad_weight = volume_quadrature.weight(iquad);
524 dealii::Point<dim,real2> phys_coord;
525 for (
int d=0;d<dim;++d) { phys_coord[d] = 0.0;}
526 std::array< dealii::Tensor<1,dim,real2>, dim > coord_grad;
527 dealii::Tensor<2,dim,real2> metric_jacobian;
528 for (
unsigned int idof=0; idof<n_metric_dofs_cell; ++idof) {
529 const unsigned int axis = fe_metric.system_to_component_index(idof).first;
530 phys_coord[axis] += coords_coeff[idof] * fe_metric.shape_value(idof, ref_point);
531 coord_grad[axis] += coords_coeff[idof] * fe_metric.shape_grad (idof, ref_point);
533 for (
int row=0;row<dim;++row) {
534 for (
int col=0;col<dim;++col) {
535 metric_jacobian[row][col] = coord_grad[row][col];
538 const real2 jacobian_determinant = dealii::determinant(metric_jacobian);
539 dealii::Tensor<2,dim,real2> jacobian_transpose_inverse;
540 jacobian_transpose_inverse = dealii::transpose(dealii::invert(metric_jacobian));
543 std::array<real2, nstate> soln_at_q;
545 std::array< dealii::Tensor<1,dim,real2>, nstate > soln_grad_at_q;
546 for (
unsigned int idof=0; idof<n_soln_dofs_cell; ++idof) {
547 const unsigned int istate = fe_solution.system_to_component_index(idof).first;
548 if (uses_solution_values) {
549 soln_at_q[istate] += soln_coeff[idof] * fe_solution.shape_value(idof,ref_point);
551 if (uses_solution_gradient) {
552 const dealii::Tensor<1,dim,real2> phys_shape_grad = dealii::contract<1,0>(jacobian_transpose_inverse, fe_solution.shape_grad(idof,ref_point));
553 soln_grad_at_q[istate] += soln_coeff[idof] * phys_shape_grad;
556 real2 volume_integrand = this->evaluate_volume_integrand(physics, phys_coord, soln_at_q, soln_grad_at_q);
558 volume_local_sum += volume_integrand * jacobian_determinant * quad_weight;
560 return volume_local_sum;
563 template <
int dim,
int nstate,
typename real,
typename MeshType>
564 template <
typename real2>
567 const unsigned int boundary_id,
568 const std::vector< real2 > &soln_coeff,
569 const dealii::FESystem<dim> &fe_solution,
570 const std::vector< real2 > &coords_coeff,
571 const dealii::FESystem<dim> &fe_metric,
572 const unsigned int face_number,
573 const dealii::Quadrature<dim-1> &fquadrature)
const 575 const unsigned int n_face_quad_pts = fquadrature.size();
576 const unsigned int n_soln_dofs_cell = soln_coeff.size();
577 const unsigned int n_metric_dofs_cell = coords_coeff.size();
579 const dealii::Quadrature<dim> face_quadrature = dealii::QProjector<dim>::project_to_face( dealii::ReferenceCell::get_hypercube(dim),
582 const dealii::Tensor<1,dim,real> surface_unit_normal = dealii::GeometryInfo<dim>::unit_normal_vector[face_number];
584 real2 boundary_local_sum = 0.0;
585 for (
unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
587 const dealii::Point<dim,double> &ref_point = face_quadrature.point(iquad);
588 const double quad_weight = face_quadrature.weight(iquad);
592 dealii::Point<dim,real2> phys_coord;
593 for (
int d=0;d<dim;++d) { phys_coord[d] = 0.0;}
594 std::array< dealii::Tensor<1,dim,real2>, dim > coord_grad;
595 dealii::Tensor<2,dim,real2> metric_jacobian;
596 for (
unsigned int idof=0; idof<n_metric_dofs_cell; ++idof) {
597 const unsigned int axis = fe_metric.system_to_component_index(idof).first;
598 phys_coord[axis] += coords_coeff[idof] * fe_metric.shape_value(idof, ref_point);
599 coord_grad[axis] += coords_coeff[idof] * fe_metric.shape_grad (idof, ref_point);
601 for (
int row=0;row<dim;++row) {
602 for (
int col=0;col<dim;++col) {
603 metric_jacobian[row][col] = coord_grad[row][col];
606 const real2 jacobian_determinant = dealii::determinant(metric_jacobian);
607 dealii::Tensor<2,dim,real2> jacobian_transpose_inverse;
608 jacobian_transpose_inverse = dealii::transpose(dealii::invert(metric_jacobian));
610 const dealii::Tensor<1,dim,real2> phys_normal = vmult(jacobian_transpose_inverse, surface_unit_normal);
611 const real2 area = norm(phys_normal);
612 const dealii::Tensor<1,dim,real2> phys_unit_normal = phys_normal/area;
614 real2 surface_jacobian_determinant = area*jacobian_determinant;
617 std::array<real2, nstate> soln_at_q;
619 std::array< dealii::Tensor<1,dim,real2>, nstate > soln_grad_at_q;
620 for (
unsigned int idof=0; idof<n_soln_dofs_cell; ++idof) {
621 const unsigned int istate = fe_solution.system_to_component_index(idof).first;
622 if (uses_solution_values) {
623 soln_at_q[istate] += soln_coeff[idof] * fe_solution.shape_value(idof,ref_point);
625 if (uses_solution_gradient) {
626 const dealii::Tensor<1,dim,real2> phys_shape_grad = dealii::contract<1,0>(jacobian_transpose_inverse, fe_solution.shape_grad(idof,ref_point));
627 soln_grad_at_q[istate] += soln_coeff[idof] * phys_shape_grad;
630 real2 boundary_integrand = this->evaluate_boundary_integrand(physics, boundary_id, phys_coord, phys_unit_normal, soln_at_q, soln_grad_at_q);
632 boundary_local_sum += boundary_integrand * surface_jacobian_determinant * quad_weight;
634 return boundary_local_sum;
637 template <
int dim,
int nstate,
typename real,
typename MeshType>
640 const unsigned int boundary_id,
641 const std::vector< real > &soln_coeff,
642 const dealii::FESystem<dim> &fe_solution,
643 const std::vector< real > &coords_coeff,
644 const dealii::FESystem<dim> &fe_metric,
645 const unsigned int face_number,
646 const dealii::Quadrature<dim-1> &fquadrature)
const 648 return evaluate_boundary_cell_functional<real>(physics, boundary_id, soln_coeff, fe_solution, coords_coeff, fe_metric, face_number, fquadrature);
651 template <
int dim,
int nstate,
typename real,
typename MeshType>
654 const unsigned int boundary_id,
655 const std::vector< Sacado::Fad::DFad<Sacado::Fad::DFad<real>> > &soln_coeff,
656 const dealii::FESystem<dim> &fe_solution,
657 const std::vector< Sacado::Fad::DFad<Sacado::Fad::DFad<real>> > &coords_coeff,
658 const dealii::FESystem<dim> &fe_metric,
659 const unsigned int face_number,
660 const dealii::Quadrature<dim-1> &fquadrature)
const 662 return evaluate_boundary_cell_functional<Sacado::Fad::DFad<Sacado::Fad::DFad<real>>>(physics, boundary_id, soln_coeff, fe_solution, coords_coeff, fe_metric, face_number, fquadrature);
665 template <
int dim,
int nstate,
typename real,
typename MeshType>
668 const std::vector< real > &soln_coeff,
669 const dealii::FESystem<dim> &fe_solution,
670 const std::vector< real > &coords_coeff,
671 const dealii::FESystem<dim> &fe_metric,
672 const dealii::Quadrature<dim> &volume_quadrature)
const 674 return evaluate_volume_cell_functional<real>(physics, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
677 template <
int dim,
int nstate,
typename real,
typename MeshType>
679 const Physics::PhysicsBase<dim,nstate,Sacado::Fad::DFad<Sacado::Fad::DFad<real>>> &physics_fad_fad,
680 const std::vector< Sacado::Fad::DFad<Sacado::Fad::DFad<real>> > &soln_coeff,
681 const dealii::FESystem<dim> &fe_solution,
682 const std::vector< Sacado::Fad::DFad<Sacado::Fad::DFad<real>> > &coords_coeff,
683 const dealii::FESystem<dim> &fe_metric,
684 const dealii::Quadrature<dim> &volume_quadrature)
const 686 return evaluate_volume_cell_functional<Sacado::Fad::DFad<Sacado::Fad::DFad<real>>>(physics_fad_fad, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
689 template <
int dim,
int nstate,
typename real,
typename MeshType>
693 pcout <<
" with value...";
695 if (dg->solution.size() == solution_value.size()
696 && dg->high_order_grid->volume_nodes.size() == volume_nodes_value.size()) {
698 auto diff_sol = dg->solution;
699 diff_sol -= solution_value;
700 const double l2_norm_sol = diff_sol.l2_norm();
702 if (l2_norm_sol == 0.0) {
704 auto diff_node = dg->high_order_grid->volume_nodes;
705 diff_node -= volume_nodes_value;
706 const double l2_norm_node = diff_node.l2_norm();
708 if (l2_norm_node == 0.0) {
709 pcout <<
" which is already assembled...";
710 compute_value =
false;
714 solution_value = dg->solution;
715 volume_nodes_value = dg->high_order_grid->volume_nodes;
718 pcout <<
" with dIdW...";
720 if (dg->solution.size() == solution_dIdW.size()
721 && dg->high_order_grid->volume_nodes.size() == volume_nodes_dIdW.size()) {
723 auto diff_sol = dg->solution;
724 diff_sol -= solution_dIdW;
725 const double l2_norm_sol = diff_sol.l2_norm();
727 if (l2_norm_sol == 0.0) {
729 auto diff_node = dg->high_order_grid->volume_nodes;
730 diff_node -= volume_nodes_dIdW;
731 const double l2_norm_node = diff_node.l2_norm();
733 if (l2_norm_node == 0.0) {
734 pcout <<
" which is already assembled...";
735 compute_dIdW =
false;
739 solution_dIdW = dg->solution;
740 volume_nodes_dIdW = dg->high_order_grid->volume_nodes;
743 pcout <<
" with dIdX...";
745 if (dg->solution.size() == solution_dIdX.size()
746 && dg->high_order_grid->volume_nodes.size() == volume_nodes_dIdX.size()) {
747 auto diff_sol = dg->solution;
748 diff_sol -= solution_dIdX;
749 const double l2_norm_sol = diff_sol.l2_norm();
751 if (l2_norm_sol == 0.0) {
753 auto diff_node = dg->high_order_grid->volume_nodes;
754 diff_node -= volume_nodes_dIdX;
755 const double l2_norm_node = diff_node.l2_norm();
757 if (l2_norm_node == 0.0) {
758 pcout <<
" which is already assembled...";
759 compute_dIdX =
false;
763 solution_dIdX = dg->solution;
764 volume_nodes_dIdX = dg->high_order_grid->volume_nodes;
767 pcout <<
" with d2IdWdW, d2IdWdX, d2IdXdX...";
769 if (dg->solution.size() == solution_d2I.size()
770 && dg->high_order_grid->volume_nodes.size() == volume_nodes_d2I.size()) {
772 auto diff_sol = dg->solution;
773 diff_sol -= solution_d2I;
774 const double l2_norm_sol = diff_sol.l2_norm();
776 if (l2_norm_sol == 0.0) {
778 auto diff_node = dg->high_order_grid->volume_nodes;
779 diff_node -= volume_nodes_d2I;
780 const double l2_norm_node = diff_node.l2_norm();
782 if (l2_norm_node == 0.0) {
784 pcout <<
" which is already assembled...";
789 solution_d2I = dg->solution;
790 volume_nodes_d2I = dg->high_order_grid->volume_nodes;
794 template <
int dim,
int nstate,
typename real,
typename MeshType>
796 const bool compute_dIdW,
797 const bool compute_dIdX,
798 const bool compute_d2I)
800 using FadType = Sacado::Fad::DFad<real>;
801 using FadFadType = Sacado::Fad::DFad<FadType>;
803 bool actually_compute_value =
true;
804 bool actually_compute_dIdW = compute_dIdW;
805 bool actually_compute_dIdX = compute_dIdX;
806 bool actually_compute_d2I = compute_d2I;
808 pcout <<
"Evaluating functional... ";
809 need_compute(actually_compute_value, actually_compute_dIdW, actually_compute_dIdX, actually_compute_d2I);
812 if (!actually_compute_value && !actually_compute_dIdW && !actually_compute_dIdX && !actually_compute_d2I) {
813 return current_functional_value;
817 real local_functional = 0.0;
820 const dealii::FESystem<dim,dim> &fe_metric = dg->high_order_grid->fe_system;
821 const unsigned int n_metric_dofs_cell = fe_metric.dofs_per_cell;
822 std::vector<dealii::types::global_dof_index> cell_metric_dofs_indices(n_metric_dofs_cell);
825 const unsigned int max_dofs_per_cell = dg->dof_handler.get_fe_collection().max_dofs_per_cell();
826 std::vector<dealii::types::global_dof_index> cell_soln_dofs_indices(max_dofs_per_cell);
827 std::vector<FadFadType> soln_coeff(max_dofs_per_cell);
828 std::vector<real> local_dIdw(max_dofs_per_cell);
830 std::vector<real> local_dIdX(n_metric_dofs_cell);
832 const auto mapping = (*(dg->high_order_grid->mapping_fe_field));
833 dealii::hp::MappingCollection<dim> mapping_collection(mapping);
835 dealii::hp::FEFaceValues<dim,dim> fe_values_collection_face (mapping_collection, dg->fe_collection, dg->face_quadrature_collection, this->face_update_flags);
837 allocate_derivatives(actually_compute_dIdW, actually_compute_dIdX, actually_compute_d2I);
839 dg->solution.update_ghost_values();
840 auto metric_cell = dg->high_order_grid->dof_handler_grid.begin_active();
841 auto soln_cell = dg->dof_handler.begin_active();
842 for( ; soln_cell != dg->dof_handler.end(); ++soln_cell, ++metric_cell) {
843 if(!soln_cell->is_locally_owned())
continue;
847 const unsigned int i_fele = soln_cell->active_fe_index();
848 const unsigned int i_quad = i_fele;
851 const dealii::FESystem<dim,dim> &fe_solution = dg->fe_collection[i_fele];
852 const unsigned int n_soln_dofs_cell = fe_solution.n_dofs_per_cell();
853 cell_soln_dofs_indices.resize(n_soln_dofs_cell);
854 soln_cell->get_dof_indices(cell_soln_dofs_indices);
855 soln_coeff.resize(n_soln_dofs_cell);
858 metric_cell->get_dof_indices (cell_metric_dofs_indices);
859 std::vector< FadFadType > coords_coeff(n_metric_dofs_cell);
860 for (
unsigned int idof = 0; idof < n_metric_dofs_cell; ++idof) {
861 coords_coeff[idof] = dg->high_order_grid->volume_nodes[cell_metric_dofs_indices[idof]];
865 unsigned int n_total_indep = 0;
866 if (actually_compute_dIdW || actually_compute_d2I) n_total_indep += n_soln_dofs_cell;
867 if (actually_compute_dIdX || actually_compute_d2I) n_total_indep += n_metric_dofs_cell;
868 unsigned int i_derivative = 0;
869 for(
unsigned int idof = 0; idof < n_soln_dofs_cell; ++idof) {
870 const real val = dg->solution[cell_soln_dofs_indices[idof]];
871 soln_coeff[idof] = val;
872 if (actually_compute_dIdW || actually_compute_d2I) soln_coeff[idof].diff(i_derivative++, n_total_indep);
874 for (
unsigned int idof = 0; idof < n_metric_dofs_cell; ++idof) {
875 const real val = dg->high_order_grid->volume_nodes[cell_metric_dofs_indices[idof]];
876 coords_coeff[idof] = val;
877 if (actually_compute_dIdX || actually_compute_d2I) coords_coeff[idof].diff(i_derivative++, n_total_indep);
879 AssertDimension(i_derivative, n_total_indep);
880 if (actually_compute_d2I) {
881 unsigned int i_derivative = 0;
882 for(
unsigned int idof = 0; idof < n_soln_dofs_cell; ++idof) {
883 const real val = dg->solution[cell_soln_dofs_indices[idof]];
884 soln_coeff[idof].val() = val;
885 soln_coeff[idof].val().diff(i_derivative++, n_total_indep);
887 for (
unsigned int idof = 0; idof < n_metric_dofs_cell; ++idof) {
888 const real val = dg->high_order_grid->volume_nodes[cell_metric_dofs_indices[idof]];
889 coords_coeff[idof].val() = val;
890 coords_coeff[idof].val().diff(i_derivative++, n_total_indep);
893 AssertDimension(i_derivative, n_total_indep);
896 const dealii::Quadrature<dim> &volume_quadrature = dg->volume_quadrature_collection[i_quad];
899 FadFadType volume_local_sum = evaluate_volume_cell_functional(*physics_fad_fad, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
902 for(
unsigned int iface = 0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface){
903 auto face = soln_cell->face(iface);
905 if(face->at_boundary()){
907 const unsigned int boundary_id = face->boundary_id();
912 volume_local_sum += this->evaluate_boundary_cell_functional(*physics_fad_fad, boundary_id, soln_coeff, fe_solution, coords_coeff, fe_metric, iface, dg->face_quadrature_collection[i_quad]);
917 local_functional += volume_local_sum.val().val();
921 if (actually_compute_dIdW) {
922 local_dIdw.resize(n_soln_dofs_cell);
923 for(
unsigned int idof = 0; idof < n_soln_dofs_cell; ++idof){
924 local_dIdw[idof] = volume_local_sum.dx(i_derivative++).val();
926 dIdw.add(cell_soln_dofs_indices, local_dIdw);
928 if (actually_compute_dIdX) {
929 local_dIdX.resize(n_metric_dofs_cell);
930 for(
unsigned int idof = 0; idof < n_metric_dofs_cell; ++idof){
931 local_dIdX[idof] = volume_local_sum.dx(i_derivative++).val();
933 dIdX.add(cell_metric_dofs_indices, local_dIdX);
935 if (actually_compute_dIdW || actually_compute_dIdX) AssertDimension(i_derivative, n_total_indep);
936 if (actually_compute_d2I) {
937 std::vector<real> dWidW(n_soln_dofs_cell);
938 std::vector<real> dWidX(n_metric_dofs_cell);
939 std::vector<real> dXidX(n_metric_dofs_cell);
942 for (
unsigned int idof=0; idof<n_soln_dofs_cell; ++idof) {
943 unsigned int j_derivative = 0;
944 const FadType dWi = volume_local_sum.dx(i_derivative++);
945 for (
unsigned int jdof=0; jdof<n_soln_dofs_cell; ++jdof) {
946 dWidW[jdof] = dWi.dx(j_derivative++);
948 d2IdWdW->add(cell_soln_dofs_indices[idof], cell_soln_dofs_indices, dWidW);
950 for (
unsigned int jdof=0; jdof<n_metric_dofs_cell; ++jdof) {
951 dWidX[jdof] = dWi.dx(j_derivative++);
953 d2IdWdX->add(cell_soln_dofs_indices[idof], cell_metric_dofs_indices, dWidX);
956 for (
unsigned int idof=0; idof<n_metric_dofs_cell; ++idof) {
958 const FadType dXi = volume_local_sum.dx(i_derivative++);
960 unsigned int j_derivative = n_soln_dofs_cell;
961 for (
unsigned int jdof=0; jdof<n_metric_dofs_cell; ++jdof) {
962 dXidX[jdof] = dXi.dx(j_derivative++);
964 d2IdXdX->add(cell_metric_dofs_indices[idof], cell_metric_dofs_indices, dXidX);
967 AssertDimension(i_derivative, n_total_indep);
970 current_functional_value = dealii::Utilities::MPI::sum(local_functional, MPI_COMM_WORLD);
972 if (actually_compute_dIdW) dIdw.compress(dealii::VectorOperation::add);
973 if (actually_compute_dIdX) dIdX.compress(dealii::VectorOperation::add);
974 if (actually_compute_d2I) {
975 d2IdWdW->compress(dealii::VectorOperation::add);
976 d2IdWdX->compress(dealii::VectorOperation::add);
977 d2IdXdX->compress(dealii::VectorOperation::add);
980 return current_functional_value;
983 template <
int dim,
int nstate,
typename real,
typename MeshType>
987 const double stepsize)
990 double local_sum_old;
991 double local_sum_new;
994 dealii::LinearAlgebra::distributed::Vector<real> dIdw;
997 dealii::IndexSet locally_owned_dofs = dg.
dof_handler.locally_owned_dofs();
998 dIdw.reinit(locally_owned_dofs, MPI_COMM_WORLD);
1001 const unsigned int max_dofs_per_cell = dg.
dof_handler.get_fe_collection().max_dofs_per_cell();
1002 std::vector<dealii::types::global_dof_index> cell_soln_dofs_indices(max_dofs_per_cell);
1003 std::vector<dealii::types::global_dof_index> neighbor_dofs_indices(max_dofs_per_cell);
1004 std::vector<real> soln_coeff(max_dofs_per_cell);
1005 std::vector<real> local_dIdw(max_dofs_per_cell);
1008 dealii::hp::MappingCollection<dim> mapping_collection(mapping);
1014 auto metric_cell = dg.
high_order_grid->dof_handler_grid.begin_active();
1016 for( ; cell != dg.
dof_handler.end(); ++cell, ++metric_cell) {
1017 if(!cell->is_locally_owned())
continue;
1020 const unsigned int i_mapp = 0;
1021 const unsigned int i_fele = cell->active_fe_index();
1022 const unsigned int i_quad = i_fele;
1025 const dealii::FESystem<dim,dim> &fe_solution = dg.
fe_collection[i_fele];
1026 const unsigned int n_soln_dofs_cell = fe_solution.n_dofs_per_cell();
1027 cell_soln_dofs_indices.resize(n_soln_dofs_cell);
1028 cell->get_dof_indices(cell_soln_dofs_indices);
1029 soln_coeff.resize(n_soln_dofs_cell);
1030 for(
unsigned int idof = 0; idof < n_soln_dofs_cell; ++idof) {
1031 soln_coeff[idof] = dg.
solution[cell_soln_dofs_indices[idof]];
1035 const dealii::FESystem<dim,dim> &fe_metric = dg.
high_order_grid->fe_system;
1036 const unsigned int n_metric_dofs_cell = fe_metric.dofs_per_cell;
1037 std::vector<dealii::types::global_dof_index> cell_metric_dofs_indices(n_metric_dofs_cell);
1038 metric_cell->get_dof_indices (cell_metric_dofs_indices);
1039 std::vector<real> coords_coeff(n_metric_dofs_cell);
1040 for (
unsigned int idof = 0; idof < n_metric_dofs_cell; ++idof) {
1041 coords_coeff[idof] = dg.
high_order_grid->volume_nodes[cell_metric_dofs_indices[idof]];
1048 local_sum_old = this->evaluate_volume_cell_functional(physics, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
1051 for(
unsigned int iface = 0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface){
1052 auto face = cell->face(iface);
1054 if(face->at_boundary()){
1055 fe_values_collection_face.reinit(cell, iface, i_quad, i_mapp, i_fele);
1057 const unsigned int boundary_id = face->boundary_id();
1061 local_sum_old += this->evaluate_boundary_cell_functional(physics, boundary_id, soln_coeff, fe_solution, coords_coeff, fe_metric, iface, dg.
face_quadrature_collection[i_quad]);
1067 local_dIdw.resize(n_soln_dofs_cell);
1068 for(
unsigned int idof = 0; idof < n_soln_dofs_cell; ++idof){
1070 for(
unsigned int idof2 = 0; idof2 < n_soln_dofs_cell; ++idof2){
1071 soln_coeff[idof2] = dg.
solution[cell_soln_dofs_indices[idof2]];
1073 soln_coeff[idof] += stepsize;
1077 local_sum_new = this->evaluate_volume_cell_functional(physics, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
1080 for(
unsigned int iface = 0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface){
1081 auto face = cell->face(iface);
1083 if(face->at_boundary()){
1084 fe_values_collection_face.reinit(cell, iface, i_quad, i_mapp, i_fele);
1086 const unsigned int boundary_id = face->boundary_id();
1090 local_sum_new += this->evaluate_boundary_cell_functional(physics, boundary_id, soln_coeff, fe_solution, coords_coeff, fe_metric, iface, dg.
face_quadrature_collection[i_quad]);
1095 local_dIdw[idof] = (local_sum_new-local_sum_old)/stepsize;
1098 dIdw.add(cell_soln_dofs_indices, local_dIdw);
1101 dIdw.compress(dealii::VectorOperation::add);
1106 template <
int dim,
int nstate,
typename real,
typename MeshType>
1110 const double stepsize)
1113 double local_sum_old;
1114 double local_sum_new;
1117 dealii::LinearAlgebra::distributed::Vector<real> dIdX_FD;
1118 allocate_dIdX(dIdX_FD);
1121 const unsigned int max_dofs_per_cell = dg.
dof_handler.get_fe_collection().max_dofs_per_cell();
1122 std::vector<dealii::types::global_dof_index> cell_soln_dofs_indices(max_dofs_per_cell);
1123 std::vector<dealii::types::global_dof_index> neighbor_dofs_indices(max_dofs_per_cell);
1124 std::vector<real> soln_coeff(max_dofs_per_cell);
1125 std::vector<real> local_dIdX(max_dofs_per_cell);
1128 dealii::hp::MappingCollection<dim> mapping_collection(mapping);
1134 auto metric_cell = dg.
high_order_grid->dof_handler_grid.begin_active();
1136 for( ; cell != dg.
dof_handler.end(); ++cell, ++metric_cell) {
1137 if(!cell->is_locally_owned())
continue;
1141 const unsigned int i_fele = cell->active_fe_index();
1142 const unsigned int i_quad = i_fele;
1145 const dealii::FESystem<dim,dim> &fe_solution = dg.
fe_collection[i_fele];
1146 const unsigned int n_soln_dofs_cell = fe_solution.n_dofs_per_cell();
1147 cell_soln_dofs_indices.resize(n_soln_dofs_cell);
1148 cell->get_dof_indices(cell_soln_dofs_indices);
1149 soln_coeff.resize(n_soln_dofs_cell);
1150 for(
unsigned int idof = 0; idof < n_soln_dofs_cell; ++idof) {
1151 soln_coeff[idof] = dg.
solution[cell_soln_dofs_indices[idof]];
1155 const dealii::FESystem<dim,dim> &fe_metric = dg.
high_order_grid->fe_system;
1156 const unsigned int n_metric_dofs_cell = fe_metric.dofs_per_cell;
1157 std::vector<dealii::types::global_dof_index> cell_metric_dofs_indices(n_metric_dofs_cell);
1158 metric_cell->get_dof_indices (cell_metric_dofs_indices);
1159 std::vector<real> coords_coeff(n_metric_dofs_cell);
1160 for (
unsigned int idof = 0; idof < n_metric_dofs_cell; ++idof) {
1161 coords_coeff[idof] = dg.
high_order_grid->volume_nodes[cell_metric_dofs_indices[idof]];
1168 local_sum_old = this->evaluate_volume_cell_functional(physics, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
1171 for(
unsigned int iface = 0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface){
1172 auto face = cell->face(iface);
1174 if(face->at_boundary()){
1176 const unsigned int boundary_id = face->boundary_id();
1181 local_sum_old += this->evaluate_boundary_cell_functional(physics, boundary_id, soln_coeff, fe_solution, coords_coeff, fe_metric, iface, dg.
face_quadrature_collection[i_quad]);
1187 local_dIdX.resize(n_metric_dofs_cell);
1188 for(
unsigned int idof = 0; idof < n_metric_dofs_cell; ++idof){
1190 for(
unsigned int idof2 = 0; idof2 < n_metric_dofs_cell; ++idof2){
1191 coords_coeff[idof2] = dg.
high_order_grid->volume_nodes[cell_metric_dofs_indices[idof2]];
1193 coords_coeff[idof] += stepsize;
1196 local_sum_new = this->evaluate_volume_cell_functional(physics, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
1199 for(
unsigned int iface = 0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface){
1200 auto face = cell->face(iface);
1202 if(face->at_boundary()){
1204 const unsigned int boundary_id = face->boundary_id();
1209 local_sum_new += this->evaluate_boundary_cell_functional(physics, boundary_id, soln_coeff, fe_solution, coords_coeff, fe_metric, iface, dg.
face_quadrature_collection[i_quad]);
1214 local_dIdX[idof] = (local_sum_new-local_sum_old)/stepsize;
1217 dIdX_FD.add(cell_metric_dofs_indices, local_dIdX);
1220 dIdX_FD.compress(dealii::VectorOperation::add);
1225 template <
int dim,
int nstate,
typename real,
typename MeshType>
1226 std::shared_ptr< Functional<dim,nstate,real,MeshType> >
1234 template <
int dim,
int nstate,
typename real,
typename MeshType>
1235 std::shared_ptr< Functional<dim,nstate,real,MeshType> >
1240 using FadFadType = Sacado::Fad::DFad<FadType>;
1246 const double normLp = param.
normLp;
1249 std::shared_ptr< ManufacturedSolutionFunction<dim,real> > weight_function_double
1251 std::shared_ptr< ManufacturedSolutionFunction<dim,FadFadType> > weight_function_adtype
1259 if(functional_type == FunctionalTypeEnum::normLp_volume){
1260 return std::make_shared<FunctionalNormLpVolume<dim,nstate,real,MeshType>>(
1265 }
else if(functional_type == FunctionalTypeEnum::normLp_boundary){
1266 return std::make_shared<FunctionalNormLpBoundary<dim,nstate,real,MeshType>>(
1273 }
else if(functional_type == FunctionalTypeEnum::weighted_integral_volume){
1274 return std::make_shared<FunctionalWeightedIntegralVolume<dim,nstate,real,MeshType>>(
1275 weight_function_double,
1276 weight_function_adtype,
1277 use_weight_function_laplacian,
1281 }
else if(functional_type == FunctionalTypeEnum::weighted_integral_boundary){
1282 return std::make_shared<FunctionalWeightedIntegralBoundary<dim,nstate,real,MeshType>>(
1283 weight_function_double,
1284 weight_function_adtype,
1285 use_weight_function_laplacian,
1291 }
else if(functional_type == FunctionalTypeEnum::error_normLp_volume){
1292 return std::make_shared<FunctionalErrorNormLpVolume<dim,nstate,real,MeshType>>(
1297 }
else if(functional_type == FunctionalTypeEnum::error_normLp_boundary){
1298 return std::make_shared<FunctionalErrorNormLpBoundary<dim,nstate,real,MeshType>>(
1305 }
else if(functional_type == FunctionalTypeEnum::lift){
1306 if constexpr(dim==2 &&
1308 std::is_same<MeshType, dealii::parallel::distributed::Triangulation<dim>>::value)
1310 return std::make_shared<LiftDragFunctional<dim,nstate,real,MeshType>>(
1314 }
else if(functional_type == FunctionalTypeEnum::drag){
1315 if constexpr(dim==2 &&
1317 std::is_same<MeshType, dealii::parallel::distributed::Triangulation<dim>>::value)
1319 return std::make_shared<LiftDragFunctional<dim,nstate,real,MeshType>>(
1323 }
else if(functional_type == FunctionalTypeEnum::solution_integral) {
1325 return std::make_shared<SolutionIntegral<dim,nstate,real,MeshType>>(dg,dg_state->pde_physics_fad_fad,
true,
false);
1326 }
else if(functional_type == FunctionalTypeEnum::outlet_pressure_integral) {
1327 if constexpr (dim==2 && nstate==dim+2){
1328 return std::make_shared<OutletPressureIntegral<dim,nstate,real,MeshType>>(dg,
true,
false);
1331 std::cout <<
"Invalid Functional." << std::endl;
ManufacturedSolutionType
Selects the manufactured solution to be used if use_manufactured_source_term=true.
Weighted volume norm functional class.
FunctionalType functional_type
Selection of functinal type.
Lp volume norm functional class.
Sacado::Fad::DFad< FadType > FadFadType
Sacado AD type that allows 2nd derivatives.
Sacado::Fad::DFad< real > FadType
Sacado AD type for first derivatives.
Base class from which Advection, Diffusion, ConvectionDiffusion, and Euler is derived.
FunctionalType
Choices for functional types to be used.
ManufacturedSolutionEnum weight_function_type
Choice of manufactured solution function to be used in weighting expression.
std::vector< unsigned int > boundary_vector
Boundary of vector ids to be considered for boundary functional evaluation.
Sacado::Fad::DFad< FadType > FadFadType
Sacado AD type that allows 2nd derivatives.
std::shared_ptr< Physics::PhysicsBase< dim, nstate, Sacado::Fad::DFad< real > > > physics
Problem physics (for calling the functional class)
Files for the baseline physics.
bool use_all_boundaries
Flag for use of all domain boundaries.
GridRefinementStudyParam grid_refinement_study_param
Contains the parameters for grid refinement study.
Manufactured solution function factory.
Functional to take the integral of the solution.
Main parameter class that contains the various other sub-parameter classes.
Lp boundary norm functional class.
dealii::LinearAlgebra::distributed::Vector< double > solution
Current modal coefficients of the solution.
Factory class to construct default functional types.
dealii::hp::QCollection< dim-1 > face_quadrature_collection
Quadrature used to evaluate face integrals.
double normLp
Choice of Lp norm exponent used in functional calculation.
FunctionalParam functional_param
Functional parameters to be used with grid refinement study.
static std::shared_ptr< ModelBase< dim, nstate, real > > create_Model(const Parameters::AllParameters *const parameters_input)
Factory to return the correct model given input parameters.
Euler equations. Derived from PhysicsBase.
std::shared_ptr< HighOrderGrid< dim, real, MeshType > > high_order_grid
High order grid that will provide the MappingFEField.
bool use_weight_function_laplacian
Flag to use weight function laplacian.
Parameterse related to the functional object.
Lp volume error norm functional class.
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.
Abstract class templated on the number of state variables.
Functional(std::shared_ptr< PHiLiP::DGBase< dim, real, MeshType >> _dg, const bool _uses_solution_values=true, const bool _uses_solution_gradient=true)
real2 compute_pressure(const std::array< real2, nstate > &conservative_soln) const
Evaluate pressure from conservative variables.
std::shared_ptr< DGBase< dim, real, MeshType > > dg
DG class pointer.
dealii::DoFHandler< dim > dof_handler
Finite Element Collection to represent the high-order grid.
Weighted boundary norm functional class.
real2 evaluate_boundary_integrand(const PHiLiP::Physics::PhysicsBase< dim, nstate, real2 > &physics, const unsigned int boundary_id, const dealii::Point< dim, real2 > &phys_coord, const dealii::Tensor< 1, dim, real2 > &normal, const std::array< real2, nstate > &soln_at_q, const std::array< dealii::Tensor< 1, dim, real2 >, nstate > &soln_grad_at_q) const
Templated boundary integrand.
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.
DGBase is independent of the number of state variables.
dealii::hp::QCollection< dim > volume_quadrature_collection
Finite Element Collection to represent the high-order grid.
Lp boundary error norm functional class.
std::shared_ptr< ManufacturedSolutionFunction< dim, real > > manufactured_solution_function
Manufactured solution function.
const dealii::hp::FECollection< dim > fe_collection
Finite Element Collection for p-finite-element to represent the solution.
dealii::ConditionalOStream pcout
Parallel std::cout that only outputs on mpi_rank==0.