[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
functional.cpp
1 // includes
2 #include "functional.h"
3 
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>
13 
14 #include <Sacado.hpp>
15 #include <algorithm>
16 #include <vector>
17 
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"
24 
26 
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)
32 {
33  dealii::Tensor<1,dim,real1> y;
34  for (int row=0;row<dim;++row) {
35  y[row] = 0.0;
36  for (int col=0;col<dim;++col) {
37  y[row] += A[row][col] * x[col];
38  }
39  }
40  return y;
41 }
42 
44 
48 template<int dim, typename real1>
49 real1 norm(const dealii::Tensor<1,dim,real1> x)
50 {
51  real1 val = 0.0;
52  for (int row=0;row<dim;++row) {
53  val += x[row] * x[row];
54  }
55  return sqrt(val);
56 }
57 
58 namespace PHiLiP {
59 
60 template <int dim, int nstate, typename real, typename MeshType>
61 FunctionalNormLpVolume<dim,nstate,real,MeshType>::FunctionalNormLpVolume(
62  const double _normLp,
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),
67  normLp(_normLp) {}
68 
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> & /*phys_coord*/,
74  const std::array<real2,nstate> & soln_at_q,
75  const std::array<dealii::Tensor<1,dim,real2>,nstate> &/*soln_grad_at_q*/) const
76 {
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);
80  return lpnorm_value;
81 }
82 
83 template <int dim, int nstate, typename real, typename MeshType>
84 FunctionalNormLpBoundary<dim,nstate,real,MeshType>::FunctionalNormLpBoundary(
85  const double _normLp,
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),
92  normLp(_normLp),
93  boundary_vector(_boundary_vector),
94  use_all_boundaries(_use_all_boundaries) {}
95 
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> & /* phys_coord */,
102  const dealii::Tensor<1,dim,real2> & /* normal */,
103  const std::array<real2,nstate> & soln_at_q,
104  const std::array<dealii::Tensor<1,dim,real2>,nstate> &/* soln_grad_at_q */) const
105 {
106  real2 lpnorm_value = 0;
107 
108  // condition for whether the current cell should be evaluated
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();
111 
112  if(!eval_boundary)
113  return lpnorm_value;
114 
115  for(unsigned int istate = 0; istate < nstate; ++istate)
116  lpnorm_value += pow(abs(soln_at_q[istate]), this->normLp);
117 
118  return lpnorm_value;
119 }
120 
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) {}
133 
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> & /*soln_grad_at_q*/,
141  std::shared_ptr<ManufacturedSolutionFunction<dim,real2>> weight_function) const
142 {
143  real2 val = 0;
144 
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));
148  }else{
149  for(unsigned int istate = 0; istate < nstate; ++istate)
150  val += soln_at_q[istate] * weight_function->value(phys_coord, istate);
151  }
152 
153  return val;
154 }
155 
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) {}
172 
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> & /* normal */,
180  const std::array<real2,nstate> & soln_at_q,
181  const std::array<dealii::Tensor<1,dim,real2>,nstate> &/* soln_grad_at_q */,
182  std::shared_ptr<ManufacturedSolutionFunction<dim,real2>> weight_function) const
183 {
184  real2 val = 0;
185 
186  // condition for whether the current cell should be evaluated
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();
189 
190  if(!eval_boundary)
191  return val;
192 
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));
196  }else{
197  for(unsigned int istate = 0; istate < nstate; ++istate)
198  val += soln_at_q[istate] * weight_function->value(phys_coord, istate);
199  }
200 
201  return val;
202 }
203 
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),
211  normLp(_normLp) {}
212 
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> &/*soln_grad_at_q*/) const
220 {
221  real2 lpnorm_value = 0;
222  for(unsigned int istate = 0; istate < nstate; ++istate){
223  const real2 uexact = physics.manufactured_solution_function->value(phys_coord, istate);
224  lpnorm_value += pow(abs(soln_at_q[istate] - uexact), this->normLp);
225  }
226  return lpnorm_value;
227 }
228 
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),
238  normLp(_normLp),
239  boundary_vector(_boundary_vector),
240  use_all_boundaries(_use_all_boundaries) {}
241 
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> & /* normal */,
249  const std::array<real2,nstate> & soln_at_q,
250  const std::array<dealii::Tensor<1,dim,real2>,nstate> &/* soln_grad_at_q */) const
251 {
252  real2 lpnorm_value = 0;
253 
254  // condition for whether the current cell should be evaluated
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();
257 
258  if(!eval_boundary)
259  return lpnorm_value;
260 
261  for(int istate = 0; istate < nstate; ++istate){
262  const real2 uexact = physics.manufactured_solution_function->value(phys_coord, istate);
263  lpnorm_value += pow(abs(soln_at_q[istate] - uexact), this->normLp);
264  }
265 
266  return lpnorm_value;
267 }
268 
269 
270 template <int dim, int nstate, typename real, typename MeshType>
271 template <typename real2>
274  const dealii::Point<dim,real2> &/*phys_coord*/,
275  const std::array<real2,nstate> &soln_at_q,
276  const std::array<dealii::Tensor<1,dim,real2>,nstate> &/*soln_grad_at_q*/) const
277 {
278  real2 val = 0;
279 
280  // integrating over the domain
281  for (int istate=0; istate<nstate; ++istate) {
282  val += soln_at_q[istate];
283  }
284 
285  return val;
286 }
287 
288 template <int dim, int nstate, typename real, typename MeshType>
290  std::shared_ptr<PHiLiP::DGBase<dim,real, MeshType>> dg_input,
291  const bool uses_solution_values,
292  const bool uses_solution_gradient)
293  : Functional<dim,nstate,real,MeshType>::Functional(dg_input,uses_solution_values,uses_solution_gradient)
294 {}
295 
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> & /*phys_coord*/,
302  const dealii::Tensor<1,dim,real2> & /*normal*/ ,
303  const std::array<real2,nstate> & soln_at_q,
304  const std::array<dealii::Tensor<1,dim,real2>,nstate> & /*soln_grad_at_q*/ ) const
305 {
306  real2 pressure = 0;
307 
308  if(boundary_id == 1002){
309  // casting it to a physics euler as it is needed for the pressure computation
310  const Physics::Euler<dim,nstate,real2>& euler_physics = dynamic_cast<const Physics::Euler<dim,nstate,real2>&>(physics);
311 
312  // converting the state vector to the point pressure
313  pressure = euler_physics.compute_pressure(soln_at_q);
314  }
315 
316  return pressure;
317 }
318 
319 
320 template <int dim, int nstate, typename real, typename MeshType>
322  std::shared_ptr<DGBase<dim,real,MeshType>> _dg,
323  const bool _uses_solution_values,
324  const bool _uses_solution_gradient)
325  : dg(_dg)
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)
332 {
333  using FadType = Sacado::Fad::DFad<real>;
334  using FadFadType = Sacado::Fad::DFad<FadType>;
335  std::shared_ptr<Physics::ModelBase<dim,nstate,FadFadType>> model_fad_fad = Physics::ModelFactory<dim,nstate,FadFadType>::create_Model(dg->all_parameters);
336  physics_fad_fad = Physics::PhysicsFactory<dim,nstate,FadFadType>::create_Physics(dg->all_parameters,model_fad_fad);
337 
338  init_vectors();
339 }
340 template <int dim, int nstate, typename real, typename MeshType>
342 {
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;
347 
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;
352 
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;
357 
358  solution_d2I.reinit(dg->solution);
359  solution_d2I *= 0.0;
360  volume_nodes_d2I.reinit(dg->high_order_grid->volume_nodes);
361  volume_nodes_d2I *= 0.0;
362 }
363 
364 template <int dim, int nstate, typename real, typename MeshType>
366  std::shared_ptr<PHiLiP::DGBase<dim,real,MeshType>> _dg,
367  std::shared_ptr<PHiLiP::Physics::PhysicsBase<dim,nstate,Sacado::Fad::DFad<Sacado::Fad::DFad<real>> >> _physics_fad_fad,
368  const bool _uses_solution_values,
369  const bool _uses_solution_gradient)
370  : Functional(_dg, _uses_solution_values, _uses_solution_gradient)
371 {
372  physics_fad_fad = _physics_fad_fad;
373 }
374 
375 template <int dim, int nstate, typename real, typename MeshType>
376 void Functional<dim,nstate,real,MeshType>::set_state(const dealii::LinearAlgebra::distributed::Vector<real> &solution_set)
377 {
378  dg->solution = solution_set;
379 }
380 
381 template <int dim, int nstate, typename real, typename MeshType>
382 void Functional<dim,nstate,real,MeshType>::set_geom(const dealii::LinearAlgebra::distributed::Vector<real> &volume_nodes_set)
383 {
384  dg->high_order_grid->volume_nodes = volume_nodes_set;
385 }
386 
387 template <int dim, int nstate, typename real, typename MeshType>
388 void Functional<dim,nstate,real,MeshType>::allocate_dIdX(dealii::LinearAlgebra::distributed::Vector<real> &dIdX) const
389 {
390  // allocating the vector
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);
397 }
398 
399 template <int dim, int nstate, typename real, typename MeshType>
400 void Functional<dim,nstate,real,MeshType>::allocate_derivatives(const bool compute_dIdW, const bool compute_dIdX, const bool compute_d2I)
401 {
402  if (compute_dIdW) {
403  // allocating the vector
404  dealii::IndexSet locally_owned_dofs = dg->dof_handler.locally_owned_dofs();
405  dIdw.reinit(locally_owned_dofs, MPI_COMM_WORLD);
406  }
407  if (compute_dIdX) {
408  allocate_dIdX(dIdX);
409  }
410  if (compute_d2I) {
411  {
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);
416  }
417 
418  {
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);
423  }
424 
425  {
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);
430  }
431  }
432 }
433 
434 
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)
441 {
442  using FadType = Sacado::Fad::DFad<real>;
443 
444  const unsigned int n_total_indep = volume_local_sum.size();
445  (void) n_total_indep; // Not used apart from assert.
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;
449 
450  if (compute_dIdW) {
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();
454  }
455  dIdw.add(cell_soln_dofs_indices, local_dIdw);
456  }
457  if (compute_dIdX) {
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();
461  }
462  dIdX.add(cell_metric_dofs_indices, local_dIdX);
463  }
464  if (compute_dIdW || compute_dIdX) AssertDimension(i_derivative, n_total_indep);
465  if (compute_d2I) {
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);
469 
470 
471  i_derivative = 0;
472  for (unsigned int idof=0; idof<n_soln_dofs_cell; ++idof) {
473 
474  unsigned int j_derivative = 0;
475  const FadType dWi = volume_local_sum.dx(i_derivative++);
476 
477  for (unsigned int jdof=0; jdof<n_soln_dofs_cell; ++jdof) {
478  dWidW[jdof] = dWi.dx(j_derivative++);
479  }
480  d2IdWdW->add(cell_soln_dofs_indices[idof], cell_soln_dofs_indices, dWidW);
481 
482  for (unsigned int jdof=0; jdof<n_metric_dofs_cell; ++jdof) {
483  dWidX[jdof] = dWi.dx(j_derivative++);
484  }
485  d2IdWdX->add(cell_soln_dofs_indices[idof], cell_metric_dofs_indices, dWidX);
486  }
487 
488  for (unsigned int idof=0; idof<n_metric_dofs_cell; ++idof) {
489 
490  const FadType dXi = volume_local_sum.dx(i_derivative++);
491 
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++);
495  }
496  d2IdXdX->add(cell_metric_dofs_indices[idof], cell_metric_dofs_indices, dXidX);
497  }
498  }
499  AssertDimension(i_derivative, n_total_indep);
500 }
501 
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
511 {
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();
515 
516  real2 volume_local_sum = 0.0;
517  for (unsigned int iquad=0; iquad<n_vol_quad_pts; ++iquad) {
518 
519  const dealii::Point<dim,double> &ref_point = volume_quadrature.point(iquad);
520  const double quad_weight = volume_quadrature.weight(iquad);
521 
522  // Obtain physical quadrature coordinates (Might be used if there is a source term or a wall distance)
523  // and evaluate metric terms such as the metric Jacobian, its inverse transpose, and its determinant
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; // Tensor initialize with zeros
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);
532  }
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];
536  }
537  }
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));
541 
542  // Evaluate the solution and gradient at the quadrature points
543  std::array<real2, nstate> soln_at_q;
544  soln_at_q.fill(0.0);
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);
550  }
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;
554  }
555  }
556  real2 volume_integrand = this->evaluate_volume_integrand(physics, phys_coord, soln_at_q, soln_grad_at_q);
557 
558  volume_local_sum += volume_integrand * jacobian_determinant * quad_weight;
559  }
560  return volume_local_sum;
561 }
562 
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
574 {
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();
578 
579  const dealii::Quadrature<dim> face_quadrature = dealii::QProjector<dim>::project_to_face( dealii::ReferenceCell::get_hypercube(dim),
580  fquadrature,
581  face_number);
582  const dealii::Tensor<1,dim,real> surface_unit_normal = dealii::GeometryInfo<dim>::unit_normal_vector[face_number];
583 
584  real2 boundary_local_sum = 0.0;
585  for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
586 
587  const dealii::Point<dim,double> &ref_point = face_quadrature.point(iquad);
588  const double quad_weight = face_quadrature.weight(iquad);
589 
590  // Obtain physical quadrature coordinates (Might be used if there is a source term or a wall distance)
591  // and evaluate metric terms such as the metric Jacobian, its inverse transpose, and its determinant
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; // Tensor initialize with zeros
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);
600  }
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];
604  }
605  }
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));
609 
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;
613 
614  real2 surface_jacobian_determinant = area*jacobian_determinant;
615 
616  // Evaluate the solution and gradient at the quadrature points
617  std::array<real2, nstate> soln_at_q;
618  soln_at_q.fill(0.0);
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);
624  }
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;
628  }
629  }
630  real2 boundary_integrand = this->evaluate_boundary_integrand(physics, boundary_id, phys_coord, phys_unit_normal, soln_at_q, soln_grad_at_q);
631 
632  boundary_local_sum += boundary_integrand * surface_jacobian_determinant * quad_weight;
633  }
634  return boundary_local_sum;
635 }
636 
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
647 {
648  return evaluate_boundary_cell_functional<real>(physics, boundary_id, soln_coeff, fe_solution, coords_coeff, fe_metric, face_number, fquadrature);
649 }
650 
651 template <int dim, int nstate, typename real, typename MeshType>
653  const Physics::PhysicsBase<dim,nstate,Sacado::Fad::DFad<Sacado::Fad::DFad<real>>> &physics,
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
661 {
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);
663 }
664 
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
673 {
674  return evaluate_volume_cell_functional<real>(physics, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
675 }
676 
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
685 {
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);
687 }
688 
689 template <int dim, int nstate, typename real, typename MeshType>
690 void Functional<dim,nstate,real,MeshType>::need_compute(bool &compute_value, bool &compute_dIdW, bool &compute_dIdX, bool &compute_d2I)
691 {
692  if (compute_value) {
693  pcout << " with value...";
694 
695  if (dg->solution.size() == solution_value.size()
696  && dg->high_order_grid->volume_nodes.size() == volume_nodes_value.size()) {
697 
698  auto diff_sol = dg->solution;
699  diff_sol -= solution_value;
700  const double l2_norm_sol = diff_sol.l2_norm();
701 
702  if (l2_norm_sol == 0.0) {
703 
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();
707 
708  if (l2_norm_node == 0.0) {
709  pcout << " which is already assembled...";
710  compute_value = false;
711  }
712  }
713  }
714  solution_value = dg->solution;
715  volume_nodes_value = dg->high_order_grid->volume_nodes;
716  }
717  if (compute_dIdW) {
718  pcout << " with dIdW...";
719 
720  if (dg->solution.size() == solution_dIdW.size()
721  && dg->high_order_grid->volume_nodes.size() == volume_nodes_dIdW.size()) {
722 
723  auto diff_sol = dg->solution;
724  diff_sol -= solution_dIdW;
725  const double l2_norm_sol = diff_sol.l2_norm();
726 
727  if (l2_norm_sol == 0.0) {
728 
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();
732 
733  if (l2_norm_node == 0.0) {
734  pcout << " which is already assembled...";
735  compute_dIdW = false;
736  }
737  }
738  }
739  solution_dIdW = dg->solution;
740  volume_nodes_dIdW = dg->high_order_grid->volume_nodes;
741  }
742  if (compute_dIdX) {
743  pcout << " with dIdX...";
744 
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();
750 
751  if (l2_norm_sol == 0.0) {
752 
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();
756 
757  if (l2_norm_node == 0.0) {
758  pcout << " which is already assembled...";
759  compute_dIdX = false;
760  }
761  }
762  }
763  solution_dIdX = dg->solution;
764  volume_nodes_dIdX = dg->high_order_grid->volume_nodes;
765  }
766  if (compute_d2I) {
767  pcout << " with d2IdWdW, d2IdWdX, d2IdXdX...";
768 
769  if (dg->solution.size() == solution_d2I.size()
770  && dg->high_order_grid->volume_nodes.size() == volume_nodes_d2I.size()) {
771 
772  auto diff_sol = dg->solution;
773  diff_sol -= solution_d2I;
774  const double l2_norm_sol = diff_sol.l2_norm();
775 
776  if (l2_norm_sol == 0.0) {
777 
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();
781 
782  if (l2_norm_node == 0.0) {
783 
784  pcout << " which is already assembled...";
785  compute_d2I = false;
786  }
787  }
788  }
789  solution_d2I = dg->solution;
790  volume_nodes_d2I = dg->high_order_grid->volume_nodes;
791  }
792 }
793 
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)
799 {
800  using FadType = Sacado::Fad::DFad<real>;
801  using FadFadType = Sacado::Fad::DFad<FadType>;
802 
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;
807 
808  pcout << "Evaluating functional... ";
809  need_compute(actually_compute_value, actually_compute_dIdW, actually_compute_dIdX, actually_compute_d2I);
810  pcout << std::endl;
811 
812  if (!actually_compute_value && !actually_compute_dIdW && !actually_compute_dIdX && !actually_compute_d2I) {
813  return current_functional_value;
814  }
815 
816  // Returned value
817  real local_functional = 0.0;
818 
819  // for taking the local derivatives
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);
823 
824  // setup it mostly the same as evaluating the value (with exception that local solution is also AD)
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); // for obtaining the local derivatives (to be copied back afterwards)
828  std::vector<real> local_dIdw(max_dofs_per_cell);
829 
830  std::vector<real> local_dIdX(n_metric_dofs_cell);
831 
832  const auto mapping = (*(dg->high_order_grid->mapping_fe_field));
833  dealii::hp::MappingCollection<dim> mapping_collection(mapping);
834 
835  dealii::hp::FEFaceValues<dim,dim> fe_values_collection_face (mapping_collection, dg->fe_collection, dg->face_quadrature_collection, this->face_update_flags);
836 
837  allocate_derivatives(actually_compute_dIdW, actually_compute_dIdX, actually_compute_d2I);
838 
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;
844 
845  // setting up the volume integration
846  // const unsigned int i_mapp = 0;
847  const unsigned int i_fele = soln_cell->active_fe_index();
848  const unsigned int i_quad = i_fele;
849 
850  // Get solution coefficients
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);
856 
857  // Get metric coefficients
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]];
862  }
863 
864  // Setup automatic differentiation
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);
873  }
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);
878  }
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);
886  }
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);
891  }
892  }
893  AssertDimension(i_derivative, n_total_indep);
894 
895  // Get quadrature point on reference cell
896  const dealii::Quadrature<dim> &volume_quadrature = dg->volume_quadrature_collection[i_quad];
897 
898  // Evaluate integral on the cell volume
899  FadFadType volume_local_sum = evaluate_volume_cell_functional(*physics_fad_fad, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
900 
901  // next looping over the faces of the cell checking for boundary elements
902  for(unsigned int iface = 0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface){
903  auto face = soln_cell->face(iface);
904 
905  if(face->at_boundary()){
906 
907  const unsigned int boundary_id = face->boundary_id();
908 
909  //fe_values_collection_face.reinit(soln_cell, iface, i_quad, i_mapp, i_fele);
910  //const dealii::FEFaceValues<dim,dim> &fe_values_face = fe_values_collection_face.get_present_fe_values();
911  //volume_local_sum += this->evaluate_cell_boundary(*physics_fad_fad, boundary_id, fe_values_face, soln_coeff);
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]);
913  }
914 
915  }
916 
917  local_functional += volume_local_sum.val().val();
918  // now getting the values and adding them to the derivaitve vector
919 
920  i_derivative = 0;
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();
925  }
926  dIdw.add(cell_soln_dofs_indices, local_dIdw);
927  }
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();
932  }
933  dIdX.add(cell_metric_dofs_indices, local_dIdX);
934  }
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);
940 
941  i_derivative = 0;
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++);
947  }
948  d2IdWdW->add(cell_soln_dofs_indices[idof], cell_soln_dofs_indices, dWidW);
949 
950  for (unsigned int jdof=0; jdof<n_metric_dofs_cell; ++jdof) {
951  dWidX[jdof] = dWi.dx(j_derivative++);
952  }
953  d2IdWdX->add(cell_soln_dofs_indices[idof], cell_metric_dofs_indices, dWidX);
954  }
955 
956  for (unsigned int idof=0; idof<n_metric_dofs_cell; ++idof) {
957 
958  const FadType dXi = volume_local_sum.dx(i_derivative++);
959 
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++);
963  }
964  d2IdXdX->add(cell_metric_dofs_indices[idof], cell_metric_dofs_indices, dXidX);
965  }
966  }
967  AssertDimension(i_derivative, n_total_indep);
968  }
969 
970  current_functional_value = dealii::Utilities::MPI::sum(local_functional, MPI_COMM_WORLD);
971  // compress before the return
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);
978  }
979 
980  return current_functional_value;
981 }
982 
983 template <int dim, int nstate, typename real, typename MeshType>
984 dealii::LinearAlgebra::distributed::Vector<real> Functional<dim,nstate,real,MeshType>::evaluate_dIdw_finiteDifferences(
987  const double stepsize)
988 {
989  // for taking the local derivatives
990  double local_sum_old;
991  double local_sum_new;
992 
993  // vector for storing the derivatives with respect to each DOF
994  dealii::LinearAlgebra::distributed::Vector<real> dIdw;
995 
996  // allocating the vector
997  dealii::IndexSet locally_owned_dofs = dg.dof_handler.locally_owned_dofs();
998  dIdw.reinit(locally_owned_dofs, MPI_COMM_WORLD);
999 
1000  // setup it mostly the same as evaluating the value (with exception that local solution is also AD)
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); // for obtaining the local derivatives (to be copied back afterwards)
1005  std::vector<real> local_dIdw(max_dofs_per_cell);
1006 
1007  const auto mapping = (*(dg.high_order_grid->mapping_fe_field));
1008  dealii::hp::MappingCollection<dim> mapping_collection(mapping);
1009 
1010  dealii::hp::FEValues<dim,dim> fe_values_collection_volume(mapping_collection, dg.fe_collection, dg.volume_quadrature_collection, this->volume_update_flags);
1011  dealii::hp::FEFaceValues<dim,dim> fe_values_collection_face (mapping_collection, dg.fe_collection, dg.face_quadrature_collection, this->face_update_flags);
1012 
1013  dg.solution.update_ghost_values();
1014  auto metric_cell = dg.high_order_grid->dof_handler_grid.begin_active();
1015  auto cell = dg.dof_handler.begin_active();
1016  for( ; cell != dg.dof_handler.end(); ++cell, ++metric_cell) {
1017  if(!cell->is_locally_owned()) continue;
1018 
1019  // setting up the volume integration
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;
1023 
1024  // Get solution coefficients
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]];
1032  }
1033 
1034  // Get metric coefficients
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]];
1042  }
1043 
1044  const dealii::Quadrature<dim> &volume_quadrature = dg.volume_quadrature_collection[i_quad];
1045 
1046  // adding the contribution from the current volume, also need to pass the solution vector on these points
1047  //local_sum_old = this->evaluate_volume_integrand(physics, fe_values_volume, soln_coeff);
1048  local_sum_old = this->evaluate_volume_cell_functional(physics, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
1049 
1050  // Next looping over the faces of the cell checking for boundary elements
1051  for(unsigned int iface = 0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface){
1052  auto face = cell->face(iface);
1053 
1054  if(face->at_boundary()){
1055  fe_values_collection_face.reinit(cell, iface, i_quad, i_mapp, i_fele);
1056 
1057  const unsigned int boundary_id = face->boundary_id();
1058 
1059  //const dealii::FEFaceValues<dim,dim> &fe_values_face = fe_values_collection_face.get_present_fe_values();
1060  //local_sum_old += this->evaluate_cell_boundary(physics, boundary_id, fe_values_face, soln_coeff);
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]);
1062  }
1063 
1064  }
1065 
1066  // now looping over all the DOFs in this cell and taking the FD
1067  local_dIdw.resize(n_soln_dofs_cell);
1068  for(unsigned int idof = 0; idof < n_soln_dofs_cell; ++idof){
1069  // for each dof copying the solution
1070  for(unsigned int idof2 = 0; idof2 < n_soln_dofs_cell; ++idof2){
1071  soln_coeff[idof2] = dg.solution[cell_soln_dofs_indices[idof2]];
1072  }
1073  soln_coeff[idof] += stepsize;
1074 
1075  // then peturb the idof'th value
1076  // local_sum_new = this->evaluate_volume_integrand(physics, fe_values_volume, soln_coeff);
1077  local_sum_new = this->evaluate_volume_cell_functional(physics, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
1078 
1079  // Next looping over the faces of the cell checking for boundary elements
1080  for(unsigned int iface = 0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface){
1081  auto face = cell->face(iface);
1082 
1083  if(face->at_boundary()){
1084  fe_values_collection_face.reinit(cell, iface, i_quad, i_mapp, i_fele);
1085 
1086  const unsigned int boundary_id = face->boundary_id();
1087 
1088  //const dealii::FEFaceValues<dim,dim> &fe_values_face = fe_values_collection_face.get_present_fe_values();
1089  //local_sum_new += this->evaluate_cell_boundary(physics, boundary_id, fe_values_face, soln_coeff);
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]);
1091  }
1092 
1093  }
1094 
1095  local_dIdw[idof] = (local_sum_new-local_sum_old)/stepsize;
1096  }
1097 
1098  dIdw.add(cell_soln_dofs_indices, local_dIdw);
1099  }
1100  // compress before the return
1101  dIdw.compress(dealii::VectorOperation::add);
1102 
1103  return dIdw;
1104 }
1105 
1106 template <int dim, int nstate, typename real, typename MeshType>
1107 dealii::LinearAlgebra::distributed::Vector<real> Functional<dim,nstate,real,MeshType>::evaluate_dIdX_finiteDifferences(
1110  const double stepsize)
1111 {
1112  // for taking the local derivatives
1113  double local_sum_old;
1114  double local_sum_new;
1115 
1116  // vector for storing the derivatives with respect to each DOF
1117  dealii::LinearAlgebra::distributed::Vector<real> dIdX_FD;
1118  allocate_dIdX(dIdX_FD);
1119 
1120  // setup it mostly the same as evaluating the value (with exception that local solution is also AD)
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); // for obtaining the local derivatives (to be copied back afterwards)
1125  std::vector<real> local_dIdX(max_dofs_per_cell);
1126 
1127  const auto mapping = (*(dg.high_order_grid->mapping_fe_field));
1128  dealii::hp::MappingCollection<dim> mapping_collection(mapping);
1129 
1130  dealii::hp::FEValues<dim,dim> fe_values_collection_volume(mapping_collection, dg.fe_collection, dg.volume_quadrature_collection, this->volume_update_flags);
1131  dealii::hp::FEFaceValues<dim,dim> fe_values_collection_face (mapping_collection, dg.fe_collection, dg.face_quadrature_collection, this->face_update_flags);
1132 
1133  dg.solution.update_ghost_values();
1134  auto metric_cell = dg.high_order_grid->dof_handler_grid.begin_active();
1135  auto cell = dg.dof_handler.begin_active();
1136  for( ; cell != dg.dof_handler.end(); ++cell, ++metric_cell) {
1137  if(!cell->is_locally_owned()) continue;
1138 
1139  // setting up the volume integration
1140  //const unsigned int i_mapp = 0;
1141  const unsigned int i_fele = cell->active_fe_index();
1142  const unsigned int i_quad = i_fele;
1143 
1144  // Get solution coefficients
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]];
1152  }
1153 
1154  // Get metric coefficients
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]];
1162  }
1163 
1164  const dealii::Quadrature<dim> &volume_quadrature = dg.volume_quadrature_collection[i_quad];
1165 
1166  // adding the contribution from the current volume, also need to pass the solution vector on these points
1167  //local_sum_old = this->evaluate_volume_integrand(physics, fe_values_volume, soln_coeff);
1168  local_sum_old = this->evaluate_volume_cell_functional(physics, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
1169 
1170  // Next looping over the faces of the cell checking for boundary elements
1171  for(unsigned int iface = 0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface){
1172  auto face = cell->face(iface);
1173 
1174  if(face->at_boundary()){
1175 
1176  const unsigned int boundary_id = face->boundary_id();
1177 
1178  //fe_values_collection_face.reinit(cell, iface, i_quad, i_mapp, i_fele);
1179  //const dealii::FEFaceValues<dim,dim> &fe_values_face = fe_values_collection_face.get_present_fe_values();
1180  //local_sum_old += this->evaluate_cell_boundary(physics, boundary_id, fe_values_face, soln_coeff);
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]);
1182  }
1183 
1184  }
1185 
1186  // now looping over all the DOFs in this cell and taking the FD
1187  local_dIdX.resize(n_metric_dofs_cell);
1188  for(unsigned int idof = 0; idof < n_metric_dofs_cell; ++idof){
1189  // for each dof copying the solution
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]];
1192  }
1193  coords_coeff[idof] += stepsize;
1194 
1195  // then peturb the idof'th value
1196  local_sum_new = this->evaluate_volume_cell_functional(physics, soln_coeff, fe_solution, coords_coeff, fe_metric, volume_quadrature);
1197 
1198  // Next looping over the faces of the cell checking for boundary elements
1199  for(unsigned int iface = 0; iface < dealii::GeometryInfo<dim>::faces_per_cell; ++iface){
1200  auto face = cell->face(iface);
1201 
1202  if(face->at_boundary()){
1203 
1204  const unsigned int boundary_id = face->boundary_id();
1205 
1206  //fe_values_collection_face.reinit(cell, iface, i_quad, i_mapp, i_fele);
1207  //const dealii::FEFaceValues<dim,dim> &fe_values_face = fe_values_collection_face.get_present_fe_values();
1208  //local_sum_new += this->evaluate_cell_boundary(physics, boundary_id, fe_values_face, soln_coeff);
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]);
1210  }
1211 
1212  }
1213 
1214  local_dIdX[idof] = (local_sum_new-local_sum_old)/stepsize;
1215  }
1216 
1217  dIdX_FD.add(cell_metric_dofs_indices, local_dIdX);
1218  }
1219  // compress before the return
1220  dIdX_FD.compress(dealii::VectorOperation::add);
1221 
1222  return dIdX_FD;
1223 }
1224 
1225 template <int dim, int nstate, typename real, typename MeshType>
1226 std::shared_ptr< Functional<dim,nstate,real,MeshType> >
1228  PHiLiP::Parameters::AllParameters const *const param,
1229  std::shared_ptr< PHiLiP::DGBase<dim,real,MeshType> > dg)
1230 {
1232 }
1233 
1234 template <int dim, int nstate, typename real, typename MeshType>
1235 std::shared_ptr< Functional<dim,nstate,real,MeshType> >
1238  std::shared_ptr< PHiLiP::DGBase<dim,real,MeshType> > dg)
1239 {
1240  using FadFadType = Sacado::Fad::DFad<FadType>;
1241 
1242  using FunctionalTypeEnum = Parameters::FunctionalParam::FunctionalType;
1243  using ManufacturedSolutionEnum = Parameters::ManufacturedSolutionParam::ManufacturedSolutionType;
1244  FunctionalTypeEnum functional_type = param.functional_type;
1245 
1246  const double normLp = param.normLp;
1247 
1248  ManufacturedSolutionEnum weight_function_type = param.weight_function_type;
1249  std::shared_ptr< ManufacturedSolutionFunction<dim,real> > weight_function_double
1251  std::shared_ptr< ManufacturedSolutionFunction<dim,FadFadType> > weight_function_adtype
1253 
1254  const bool use_weight_function_laplacian = param.use_weight_function_laplacian;
1255 
1256  std::vector<unsigned int> boundary_vector = param.boundary_vector;
1257  const bool use_all_boundaries = param.use_all_boundaries;
1258 
1259  if(functional_type == FunctionalTypeEnum::normLp_volume){
1260  return std::make_shared<FunctionalNormLpVolume<dim,nstate,real,MeshType>>(
1261  normLp,
1262  dg,
1263  true,
1264  false);
1265  }else if(functional_type == FunctionalTypeEnum::normLp_boundary){
1266  return std::make_shared<FunctionalNormLpBoundary<dim,nstate,real,MeshType>>(
1267  normLp,
1268  boundary_vector,
1269  use_all_boundaries,
1270  dg,
1271  true,
1272  false);
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,
1278  dg,
1279  true,
1280  false);
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,
1286  boundary_vector,
1287  use_all_boundaries,
1288  dg,
1289  true,
1290  false);
1291  }else if(functional_type == FunctionalTypeEnum::error_normLp_volume){
1292  return std::make_shared<FunctionalErrorNormLpVolume<dim,nstate,real,MeshType>>(
1293  normLp,
1294  dg,
1295  true,
1296  false);
1297  }else if(functional_type == FunctionalTypeEnum::error_normLp_boundary){
1298  return std::make_shared<FunctionalErrorNormLpBoundary<dim,nstate,real,MeshType>>(
1299  normLp,
1300  boundary_vector,
1301  use_all_boundaries,
1302  dg,
1303  true,
1304  false);
1305  }else if(functional_type == FunctionalTypeEnum::lift){
1306  if constexpr(dim==2 &&
1307  nstate==(dim+2) &&
1308  std::is_same<MeshType, dealii::parallel::distributed::Triangulation<dim>>::value)
1309  {
1310  return std::make_shared<LiftDragFunctional<dim,nstate,real,MeshType>>(
1311  dg,
1313  }
1314  }else if(functional_type == FunctionalTypeEnum::drag){
1315  if constexpr(dim==2 &&
1316  nstate==(dim+2) &&
1317  std::is_same<MeshType, dealii::parallel::distributed::Triangulation<dim>>::value)
1318  {
1319  return std::make_shared<LiftDragFunctional<dim,nstate,real,MeshType>>(
1320  dg,
1322  }
1323  }else if(functional_type == FunctionalTypeEnum::solution_integral) {
1324  std::shared_ptr< DGBaseState<dim,nstate,double,MeshType>> dg_state = std::dynamic_pointer_cast< DGBaseState<dim,nstate,double, MeshType>>(dg);
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);
1329  }
1330  }else{
1331  std::cout << "Invalid Functional." << std::endl;
1332  }
1333 
1334  return nullptr;
1335 }
1336 
1337 // dealii::Triangulation<PHILIP_DIM>
1343 
1349 
1355 
1361 
1367 
1373 
1379 
1385 
1386 // dealii::parallel::shared::Triangulation<PHILIP_DIM>
1392 
1398 
1404 
1410 
1416 
1422 
1428 
1434 
1440 
1441 #if PHILIP_DIM != 1
1442 // dealii::parallel::distributed::Triangulation<PHILIP_DIM>
1448 
1454 
1460 
1466 
1472 
1478 
1484 
1490 
1496 
1497 #endif
1498 
1499 #if PHILIP_DIM == 2
1501 #endif
1502 
1503 } // PHiLiP namespace
ManufacturedSolutionType
Selects the manufactured solution to be used if use_manufactured_source_term=true.
Weighted volume norm functional class.
Definition: functional.h:456
FunctionalType functional_type
Selection of functinal type.
Lp volume norm functional class.
Definition: functional.h:333
Sacado::Fad::DFad< FadType > FadFadType
Sacado AD type that allows 2nd derivatives.
Definition: ADTypes.hpp:12
Sacado::Fad::DFad< real > FadType
Sacado AD type for first derivatives.
Definition: functional.h:45
Base class from which Advection, Diffusion, ConvectionDiffusion, and Euler is derived.
Definition: physics.h:34
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.
Definition: functional.h:46
std::shared_ptr< Physics::PhysicsBase< dim, nstate, Sacado::Fad::DFad< real > > > physics
Problem physics (for calling the functional class)
Definition: adjoint.h:125
Files for the baseline physics.
Definition: ADTypes.hpp:10
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.
Definition: functional.h:710
Main parameter class that contains the various other sub-parameter classes.
Lp boundary norm functional class.
Definition: functional.h:389
dealii::LinearAlgebra::distributed::Vector< double > solution
Current modal coefficients of the solution.
Definition: dg_base.hpp:409
Factory class to construct default functional types.
Definition: functional.h:813
dealii::hp::QCollection< dim-1 > face_quadrature_collection
Quadrature used to evaluate face integrals.
Definition: dg_base.hpp:610
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.
Definition: euler.h:78
std::shared_ptr< HighOrderGrid< dim, real, MeshType > > high_order_grid
High order grid that will provide the MappingFEField.
Definition: dg_base.hpp:655
bool use_weight_function_laplacian
Flag to use weight function laplacian.
Parameterse related to the functional object.
Lp volume error norm functional class.
Definition: functional.h:594
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.
Definition: euler.cpp:369
std::shared_ptr< DGBase< dim, real, MeshType > > dg
DG class pointer.
Definition: adjoint.h:121
dealii::DoFHandler< dim > dof_handler
Finite Element Collection to represent the high-order grid.
Definition: dg_base.hpp:652
Weighted boundary norm functional class.
Definition: functional.h:521
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.
Definition: functional.cpp:298
Functional base class.
Definition: functional.h:43
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.
Definition: functional.cpp:272
DGBase is independent of the number of state variables.
Definition: dg_base.hpp:82
dealii::hp::QCollection< dim > volume_quadrature_collection
Finite Element Collection to represent the high-order grid.
Definition: dg_base.hpp:608
Lp boundary error norm functional class.
Definition: functional.h:651
std::shared_ptr< ManufacturedSolutionFunction< dim, real > > manufactured_solution_function
Manufactured solution function.
Definition: physics.h:71
const dealii::hp::FECollection< dim > fe_collection
Finite Element Collection for p-finite-element to represent the solution.
Definition: dg_base.hpp:597
dealii::ConditionalOStream pcout
Parallel std::cout that only outputs on mpi_rank==0.
Definition: adjoint.h:152