[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
diffusion_exact_adjoint.cpp
1 // includes
2 #include <stdlib.h>
3 #include <iostream>
4 
5 #include <deal.II/base/convergence_table.h>
6 
7 #include <deal.II/dofs/dof_tools.h>
8 
9 #include <deal.II/grid/grid_generator.h>
10 #include <deal.II/grid/grid_refinement.h>
11 #include <deal.II/grid/grid_tools.h>
12 #include <deal.II/grid/grid_out.h>
13 #include <deal.II/grid/grid_in.h>
14 
15 #include <deal.II/numerics/vector_tools.h>
16 #include <deal.II/numerics/data_out.h>
17 
18 #include <deal.II/fe/fe_values.h>
19 #include <deal.II/fe/mapping_q.h>
20 
21 #include "diffusion_exact_adjoint.h"
22 
23 #include "parameters/all_parameters.h"
24 
25 #include "physics/euler.h"
26 #include "physics/manufactured_solution.h"
27 #include "dg/dg_factory.hpp"
28 #include "dg/dg_base_state.hpp"
29 #include "ode_solver/ode_solver_factory.h"
30 
31 #include "functional/functional.h"
32 #include "functional/adjoint.h"
33 
34 #include "post_processor/physics_post_processor.h"
35 
36 #include "ADTypes.hpp"
37 
38 namespace PHiLiP {
39 namespace Tests {
40 // built own physics classes here for one time use and added a function to pass them directly to dg state
41 
42 // manufactured solution in u
43 template <int dim, typename real>
44 real ManufacturedSolutionU<dim,real>::value(const dealii::Point<dim,real> &pos, const unsigned int /*istate*/) const
45 {
46  real val = 0;
47 
48  if(dim == 1){
49  real x = pos[0];
50  val = (-1.0*pow(x,6)+3.0*pow(x,5)-3.0*pow(x,4)+pow(x,3));
51  }else if(dim == 2){
52  real x = pos[0], y = pos[1];
53  val = (-1.0*pow(x,6)+3.0*pow(x,5)-3.0*pow(x,4)+pow(x,3))
54  * (-1.0*pow(y,6)+3.0*pow(y,5)-3.0*pow(y,4)+pow(y,3));
55  }else if(dim == 3){
56  real x = pos[0], y = pos[1], z = pos[2];
57  val = (-1.0*pow(x,6)+3.0*pow(x,5)-3.0*pow(x,4)+pow(x,3))
58  * (-1.0*pow(y,6)+3.0*pow(y,5)-3.0*pow(y,4)+pow(y,3))
59  * (-1.0*pow(z,6)+3.0*pow(z,5)-3.0*pow(z,4)+pow(z,3));
60  }
61 
62  return val;
63 }
64 
65 // gradient of the solution in u
66 template <int dim, typename real>
67 dealii::Tensor<1,dim,real> ManufacturedSolutionU<dim,real>::gradient(const dealii::Point<dim,real> &pos, const unsigned int /*istate*/) const
68 {
69  dealii::Tensor<1,dim,real> gradient;
70 
71  if(dim == 1){
72  real x = pos[0];
73  gradient[0] = (-6.0*pow(x,5)+15.0*pow(x,4)-12.0*pow(x,3)+3.0*pow(x,2));
74  }else if(dim == 2){
75  real x = pos[0], y = pos[1];
76  gradient[0] = (-6.0*pow(x,5)+15.0*pow(x,4)-12.0*pow(x,3)+3.0*pow(x,2))
77  * (-1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3));
78  gradient[1] = (-1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
79  * (-6.0*pow(y,5)+15.0*pow(y,4)-12.0*pow(y,3)+3.0*pow(y,2));
80  }else if(dim == 3){
81  real x = pos[0], y = pos[1], z = pos[2];
82  gradient[0] = (-6.0*pow(x,5)+15.0*pow(x,4)-12.0*pow(x,3)+3.0*pow(x,2))
83  * (-1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
84  * (-1.0*pow(z,6)+ 3.0*pow(z,5)- 3.0*pow(z,4)+ pow(z,3));
85  gradient[1] = (-1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
86  * (-6.0*pow(y,5)+15.0*pow(y,4)-12.0*pow(y,3)+3.0*pow(y,2))
87  * (-1.0*pow(z,6)+ 3.0*pow(z,5)- 3.0*pow(z,4)+ pow(z,3));
88  gradient[2] = (-1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
89  * (-1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
90  * (-6.0*pow(z,5)+15.0*pow(z,4)-12.0*pow(z,3)+3.0*pow(z,2));
91  }
92 
93  return gradient;
94 }
95 
96 // manufactured solution in v
97 template <int dim, typename real>
98 real ManufacturedSolutionV<dim,real>::value(const dealii::Point<dim,real> &pos, const unsigned int /*istate*/) const
99 {
100  const double pi = std::acos(-1);
101 
102  real val = 0;
103 
104  if(dim == 1){
105  real x = pos[0];
106  val = (sin(pi*x));
107  }else if(dim == 2){
108  real x = pos[0], y = pos[1];
109  val = (sin(pi*x))
110  * (sin(pi*y));
111  }else if(dim == 3){
112  real x = pos[0], y = pos[1], z = pos[2];
113  val = (sin(pi*x))
114  * (sin(pi*y))
115  * (sin(pi*z));
116  }
117 
118  return val;
119 }
120 
121 // gradient of the solution in v
122 template <int dim, typename real>
123 dealii::Tensor<1,dim,real> ManufacturedSolutionV<dim,real>::gradient(const dealii::Point<dim,real> &pos, const unsigned int /*istate*/) const
124 {
125  const double pi = std::acos(-1);
126 
127  dealii::Tensor<1,dim,real> gradient;
128 
129  if(dim == 1){
130  real x = pos[0];
131  gradient[0] = pi*cos(pi*x);
132  }else if(dim == 2){
133  real x = pos[0], y = pos[1];
134  gradient[0] = (pi*cos(pi*x))
135  * ( sin(pi*y));
136  gradient[1] = ( sin(pi*x))
137  * (pi*cos(pi*y));
138  }else if(dim == 3){
139  real x = pos[0], y = pos[1], z = pos[2];
140  gradient[0] = (pi*cos(pi*x))
141  * ( sin(pi*y))
142  * ( sin(pi*z));
143  gradient[1] = ( sin(pi*x))
144  * (pi*cos(pi*y))
145  * ( sin(pi*z));
146  gradient[2] = ( sin(pi*x))
147  * ( sin(pi*y))
148  * (pi*cos(pi*z));
149  }
150 
151  return gradient;
152 }
153 
154 /* Defining the physics objects to be used */
155 template <int dim, int nstate, typename real>
157  const dealii::Point<dim,real> &pos,
158  const std::array<real,nstate> &/*solution*/,
159  const real /*current_time*/,
160  const dealii::types::global_dof_index /*cell_index*/) const
161 {
162  std::array<real,nstate> source;
163 
164  real val = 0;
165 
166  if(dim == 1){
167  real x = pos[0];
168  val = (-30.0*pow(x,4)+60.0*pow(x,3)-36.0*pow(x,2)+6.0*x);
169  }else if(dim == 2){
170  real x = pos[0], y = pos[1];
171  val = (-30.0*pow(x,4)+60.0*pow(x,3)-36.0*pow(x,2)+6.0*x)
172  * ( -1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
173  + ( -1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
174  * (-30.0*pow(y,4)+60.0*pow(y,3)-36.0*pow(y,2)+6.0*y);
175  }else if(dim == 3){
176  real x = pos[0], y = pos[1], z = pos[2];
177  val = (-30.0*pow(x,4)+60.0*pow(x,3)-36.0*pow(x,2)+6.0*x)
178  * ( -1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
179  * ( -1.0*pow(z,6)+ 3.0*pow(z,5)- 3.0*pow(z,4)+ pow(z,3))
180  + ( -1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
181  * (-30.0*pow(y,4)+60.0*pow(y,3)-36.0*pow(y,2)+6.0*y)
182  * ( -1.0*pow(z,6)+ 3.0*pow(z,5)- 3.0*pow(z,4)+ pow(z,3))
183  + ( -1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
184  * ( -1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
185  * (-30.0*pow(z,4)+60.0*pow(z,3)-36.0*pow(z,2)+6.0*z);
186  }
187 
188  for(int istate = 0; istate < nstate; ++istate)
189  source[istate] = val;
190 
191  return source;
192 }
193 
194 template <int dim, int nstate, typename real>
196  const dealii::Point<dim,real> &pos) const
197 {
198  const double pi = std::acos(-1);
199 
200  real val = 0;
201 
202  if(dim == 1){
203  real x = pos[0];
204  val = -pi*pi*sin(pi*x);
205  }else if(dim == 2){
206  real x = pos[0], y = pos[1];
207  val = -pi*pi*sin(pi*x)
208  * sin(pi*y)
209  + sin(pi*x)
210  * -pi*pi*sin(pi*y);
211  }else if(dim == 3){
212  real x = pos[0], y = pos[1], z = pos[2];
213  val = -pi*pi*sin(pi*x)
214  * sin(pi*y)
215  * sin(pi*z)
216  + sin(pi*x)
217  * -pi*pi*sin(pi*y)
218  * sin(pi*z)
219  + sin(pi*x)
220  * sin(pi*y)
221  * -pi*pi*sin(pi*z);
222  }
223 
224  return val;
225 }
226 
227 template <int dim, int nstate, typename real>
229  const dealii::Point<dim,real> &pos,
230  const std::array<real,nstate> &/*solution*/,
231  const real /*current_time*/,
232  const dealii::types::global_dof_index /*cell_index*/) const
233 {
234  const double pi = std::acos(-1);
235 
236  std::array<real,nstate> source;
237 
238  real val = 0;
239 
240  if(dim == 1){
241  real x = pos[0];
242  val = -pi*pi*sin(pi*x);
243  }else if(dim == 2){
244  real x = pos[0], y = pos[1];
245  val = -pi*pi*sin(pi*x)
246  * sin(pi*y)
247  + sin(pi*x)
248  * -pi*pi*sin(pi*y);
249  }else if(dim == 3){
250  real x = pos[0], y = pos[1], z = pos[2];
251  val = -pi*pi*sin(pi*x)
252  * sin(pi*y)
253  * sin(pi*z)
254  + sin(pi*x)
255  * -pi*pi*sin(pi*y)
256  * sin(pi*z)
257  + sin(pi*x)
258  * sin(pi*y)
259  * -pi*pi*sin(pi*z);
260  }
261 
262  for(int istate = 0; istate < nstate; ++istate)
263  source[istate] = val;
264 
265  return source;
266 }
267 
268 template <int dim, int nstate, typename real>
270  const dealii::Point<dim,real> &pos) const
271 {
272  real val = 0;
273 
274  if(dim == 1){
275  real x = pos[0];
276  val = (-30.0*pow(x,4)+60.0*pow(x,3)-36.0*pow(x,2)+6.0*x);
277  }else if(dim == 2){
278  real x = pos[0], y = pos[1];
279  val = (-30.0*pow(x,4)+60.0*pow(x,3)-36.0*pow(x,2)+6.0*x)
280  * ( -1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
281  + ( -1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
282  * (-30.0*pow(y,4)+60.0*pow(y,3)-36.0*pow(y,2)+6.0*y);
283  }else if(dim == 3){
284  real x = pos[0], y = pos[1], z = pos[2];
285  val = (-30.0*pow(x,4)+60.0*pow(x,3)-36.0*pow(x,2)+6.0*x)
286  * ( -1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
287  * ( -1.0*pow(z,6)+ 3.0*pow(z,5)- 3.0*pow(z,4)+ pow(z,3))
288  + ( -1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
289  * (-30.0*pow(y,4)+60.0*pow(y,3)-36.0*pow(y,2)+6.0*y)
290  * ( -1.0*pow(z,6)+ 3.0*pow(z,5)- 3.0*pow(z,4)+ pow(z,3))
291  + ( -1.0*pow(x,6)+ 3.0*pow(x,5)- 3.0*pow(x,4)+ pow(x,3))
292  * ( -1.0*pow(y,6)+ 3.0*pow(y,5)- 3.0*pow(y,4)+ pow(y,3))
293  * (-30.0*pow(z,4)+60.0*pow(z,3)-36.0*pow(z,2)+6.0*z);
294  }
295 
296  return val;
297 }
298 
299 // Funcitonal that performs the inner product over the entire domain
300 template <int dim, int nstate, typename real>
301 template <typename real2>
304  const dealii::Point<dim,real2> &phys_coord,
305  const std::array<real2,nstate> &soln_at_q,
306  const std::array<dealii::Tensor<1,dim,real2>,nstate> &/*soln_grad_at_q*/) const
307 {
308  real2 val = 0;
309 
310  // casting our physics object into a diffusion_objective object
311  const diffusion_objective<dim,nstate,real2> &diff_physics = dynamic_cast< const diffusion_objective<dim,nstate,real2> & >(physics);
312 
313  // evaluating the associated objective function weighting at the quadrature point
314  real2 objective_value = diff_physics.objective_function(phys_coord);
315 
316  // integrating over the domain (adding istate loop but should always be 1)
317  for (int istate=0; istate<nstate; ++istate) {
318  val += soln_at_q[istate] * objective_value;
319  }
320 
321  return val;
322 }
323 
324 
325 template <int dim, int nstate>
327  TestsBase::TestsBase(parameters_input){}
328 
329 template <int dim, int nstate>
331 {
332  pcout << "Running diffusion exact adjoint test case." << std::endl;
333 
334  // getting the problem parameters
336  using GridEnum = ManParam::GridEnum;
339 
341 
342  Assert(dim == param.dimension, dealii::ExcDimensionMismatch(dim, param.dimension));
343 
344  ManParam manu_grid_conv_param = param.manufactured_convergence_study_param;
345 
346  const unsigned int p_start = manu_grid_conv_param.degree_start;
347  const unsigned int p_end = manu_grid_conv_param.degree_end;
348  const unsigned int n_grids = manu_grid_conv_param.number_of_grids;
349 
350  // checking that the diffusion equation has been selected
351  PdeEnum pde_type = param.pde_type;
352  bool convection, diffusion;
353  if(pde_type == PdeEnum::diffusion){
354  convection = false;
355  diffusion = true;
356  }else{
357  pcout << "Can't run diffusion_exact_adjoint test case with other PDE types." << std::endl;
358  }
359 
360  // creating the physics objects
361  std::shared_ptr< Physics::PhysicsBase<dim, nstate, double> > physics_u_double
362  = std::make_shared< diffusion_u<dim, nstate, double> >(&param, convection, diffusion);
363  std::shared_ptr< Physics::PhysicsBase<dim, nstate, double> > physics_v_double
364  = std::make_shared< diffusion_v<dim, nstate, double> >(&param, convection, diffusion);
365 
366  // for adjoint, also need the AD'd physics
367  std::shared_ptr< Physics::PhysicsBase<dim, nstate, FadType> > physics_u_fadtype
368  = std::make_shared< diffusion_u<dim, nstate, FadType> >(&param, convection, diffusion);
369  std::shared_ptr< Physics::PhysicsBase<dim, nstate, FadType> > physics_v_fadtype
370  = std::make_shared< diffusion_v<dim, nstate, FadType> >(&param, convection, diffusion);
371 
372  std::shared_ptr< Physics::PhysicsBase<dim, nstate, RadType> > physics_u_radtype
373  = std::make_shared< diffusion_u<dim, nstate, RadType> >(&param, convection, diffusion);
374  std::shared_ptr< Physics::PhysicsBase<dim, nstate, RadType> > physics_v_radtype
375  = std::make_shared< diffusion_v<dim, nstate, RadType> >(&param, convection, diffusion);
376 
377  std::shared_ptr< Physics::PhysicsBase<dim, nstate, FadFadType> > physics_u_fadfadtype
378  = std::make_shared< diffusion_u<dim, nstate, FadFadType> >(&param, convection, diffusion);
379  std::shared_ptr< Physics::PhysicsBase<dim, nstate, FadFadType> > physics_v_fadfadtype
380  = std::make_shared< diffusion_v<dim, nstate, FadFadType> >(&param, convection, diffusion);
381 
382  std::shared_ptr< Physics::PhysicsBase<dim, nstate, RadFadType> > physics_u_radfadtype
383  = std::make_shared< diffusion_u<dim, nstate, RadFadType> >(&param, convection, diffusion);
384  std::shared_ptr< Physics::PhysicsBase<dim, nstate, RadFadType> > physics_v_radfadtype
385  = std::make_shared< diffusion_v<dim, nstate, RadFadType> >(&param, convection, diffusion);
386 
387  // exact value to be used in checks below
388  const double pi = std::acos(-1);
389  double exact_val = 0;
390 
391  if(dim == 1){
392  exact_val = (144*(pow(pi,2)-10)/pow(pi,5));
393  }else if(dim == 2){
394  exact_val = 2 * (-144*(pow(pi,2)-10)/pow(pi,7)) * (144*(pow(pi,2)-10)/pow(pi,5));
395  }else if(dim == 3){
396  exact_val = 3 * pow(-144*(pow(pi,2)-10)/pow(pi,7), 2) * (144*(pow(pi,2)-10)/pow(pi,5));
397  }
398 
399  // checks
400  std::vector<dealii::ConvergenceTable> convergence_table_vector;
401 
402  unsigned int n_fail = 0;
403 
404  for(unsigned int poly_degree = p_start; poly_degree <= p_end; ++poly_degree){
405  pcout << "Starting polynomial order: " << poly_degree << std::endl;
406 
407  // grid_study refines P0 additinoally here (unused in current param file)
408  std::vector<double> grid_size(n_grids);
409  std::vector<double> output_error_u(n_grids);
410  std::vector<double> output_error_v(n_grids);
411  std::vector<double> soln_error_u(n_grids);
412  std::vector<double> soln_error_v(n_grids);
413  std::vector<double> adj_error_u(n_grids);
414  std::vector<double> adj_error_v(n_grids);
415 
416  // for outputing cell-wise erorr distribution
417  dealii::Vector<double> cellError_soln_u;
418  dealii::Vector<double> cellError_soln_v;
419  dealii::Vector<double> cellError_adj_u;
420  dealii::Vector<double> cellError_adj_v;
421 
422  const std::vector<int> n_1d_cells = get_number_1d_cells(n_grids);
423 
424  dealii::ConvergenceTable convergence_table;
425 
426 #if PHILIP_DIM==1
427  using Triangulation = dealii::Triangulation<PHILIP_DIM>;
428 #else
429  using Triangulation = dealii::parallel::distributed::Triangulation<PHILIP_DIM>;
430 #endif
431  std::shared_ptr <Triangulation> grid = std::make_shared<Triangulation> (
432 #if PHILIP_DIM!=1
433  this->mpi_communicator,
434 #endif
435  typename dealii::Triangulation<dim>::MeshSmoothing(
436  dealii::Triangulation<dim>::smoothing_on_refinement |
437  dealii::Triangulation<dim>::smoothing_on_coarsening));
438 
439  // dimensions of the mesh
440  const double left = 0.0;
441  const double right = 1.0;
442 
443  for(unsigned int igrid = 0; igrid < n_grids; ++igrid){
444  // grid generation
445  grid->clear();
446  dealii::GridGenerator::subdivided_hyper_cube(*grid, n_1d_cells[igrid], left, right);
447  for (auto cell = grid->begin_active(); cell != grid->end(); ++cell) {
448  cell->set_material_id(9002);
449  for (unsigned int face=0; face<dealii::GeometryInfo<dim>::faces_per_cell; ++face) {
450  if (cell->face(face)->at_boundary()) cell->face(face)->set_boundary_id (1000);
451  }
452  }
453 
454  // since a different grid is constructed each time, need to also generate a new DG
455  // I don't think this would work outside loop since grid.clear() req's no subscriptors
456  std::shared_ptr < DGBase<dim, double> > dg_u = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree, grid);
457  std::shared_ptr < DGBase<dim, double> > dg_v = DGFactory<dim,double>::create_discontinuous_galerkin(&param, poly_degree, grid);
458 
459  // casting to dg state
460  std::shared_ptr< DGBaseState<dim,nstate,double> > dg_state_u = std::dynamic_pointer_cast< DGBaseState<dim,nstate,double> >(dg_u);
461  std::shared_ptr< DGBaseState<dim,nstate,double> > dg_state_v = std::dynamic_pointer_cast< DGBaseState<dim,nstate,double> >(dg_v);
462 
463  // now overriding the original physics on each
464  dg_state_u->set_physics(physics_u_double, physics_u_fadtype, physics_u_radtype, physics_u_fadfadtype, physics_u_radfadtype);
465  dg_state_v->set_physics(physics_v_double, physics_v_fadtype, physics_v_radtype, physics_v_fadfadtype, physics_v_radfadtype);
466 
467  dg_u->allocate_system();
468  dg_v->allocate_system();
469 
470  dg_u->solution *= 0.0;
471  dg_v->solution *= 0.0;
472  //dg_u->solution.add(1.1);
473  //dg_v->solution.add(1.1);
474 
475  // Create ODE solvers using the factory and providing the DG object
476  std::shared_ptr<ODE::ODESolverBase<dim, double>> ode_solver_u = ODE::ODESolverFactory<dim, double>::create_ODESolver(dg_u);
477  std::shared_ptr<ODE::ODESolverBase<dim, double>> ode_solver_v = ODE::ODESolverFactory<dim, double>::create_ODESolver(dg_v);
478 
479  // solving
480  ode_solver_u->steady_state();
481  ode_solver_v->steady_state();
482 
483  pcout << "Creating DiffusionFunctional... " << std::endl;
484  // functional for computations
485  auto diffusion_functional_u = std::make_shared<DiffusionFunctional<dim,nstate,double>>(dg_u,physics_u_fadfadtype,true,false);
486  auto diffusion_functional_v = std::make_shared<DiffusionFunctional<dim,nstate,double>>(dg_v,physics_v_fadfadtype,true,false);
487 
488  pcout << "Evaluating functional... " << std::endl;
489  // evaluating functionals from both methods
490  double functional_val_u = diffusion_functional_u->evaluate_functional(false,false);
491  double functional_val_v = diffusion_functional_v->evaluate_functional(false,false);
492 
493  // comparison betweent the values, add these to the convergence table
494  pcout << std::endl << "Val1 = " << functional_val_u << "\tVal2 = " << functional_val_v << std::endl << std::endl;
495 
496  // evaluating the error of this measure
497  double error_functional_u = std::abs(functional_val_u-exact_val);
498  double error_functional_v = std::abs(functional_val_v-exact_val);
499 
500  pcout << std::endl << "error_val1 = " << error_functional_u << "\terror_val2 = " << error_functional_v << std::endl << std::endl;
501 
502  // // Initializing the adjoints for each problem
503  Adjoint<dim, nstate, double> adj_u(dg_u, diffusion_functional_u, physics_u_fadtype);
504  Adjoint<dim, nstate, double> adj_v(dg_v, diffusion_functional_v, physics_v_fadtype);
505 
506  // solving for each coarse adjoint
507  pcout << "Solving for the discrete adjoints." << std::endl;
508  dealii::LinearAlgebra::distributed::Vector<double> adjoint_u = adj_u.coarse_grid_adjoint();
509  dealii::LinearAlgebra::distributed::Vector<double> adjoint_v = adj_v.coarse_grid_adjoint();
510 
511  // using overintegration for estimating the error in the adjoint
512  int overintegrate = 10;
513  dealii::QGauss<dim> quad_extra(poly_degree+overintegrate);
514  dealii::FEValues<dim,dim> fe_values_extra(*(dg_u->high_order_grid->mapping_fe_field), dg_u->fe_collection[poly_degree], quad_extra,
515  dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points);
516  const unsigned int n_quad_pts = fe_values_extra.n_quadrature_points;
517 
518  // examining the convergence of the soln compared to same manufactured solution
519  double l2error_soln_u = 0;
520  double l2error_soln_v = 0;
521 
522  // evaluate the error in the adjoints compared to the opposing manufactured solutions
523  double l2error_adj_u = 0;
524  double l2error_adj_v = 0;
525 
526  // for values at each quadrature point
527  std::array<double,nstate> soln_at_q_u;
528  std::array<double,nstate> soln_at_q_v;
529  std::array<double,nstate> adj_at_q_u;
530  std::array<double,nstate> adj_at_q_v;
531 
532  // reinit vectors for error distribution
533  cellError_soln_u.reinit(grid->n_active_cells());
534  cellError_soln_v.reinit(grid->n_active_cells());
535  cellError_adj_u.reinit(grid->n_active_cells());
536  cellError_adj_v.reinit(grid->n_active_cells());
537 
538  std::vector<dealii::types::global_dof_index> dofs_indices(fe_values_extra.dofs_per_cell);
539  for(auto cell = dg_u->dof_handler.begin_active(); cell != dg_u->dof_handler.end(); ++cell){
540  if(!cell->is_locally_owned()) continue;
541 
542  fe_values_extra.reinit (cell);
543  cell->get_dof_indices(dofs_indices);
544 
545  double cell_l2error_soln_u = 0;
546  double cell_l2error_soln_v = 0;
547  double cell_l2error_adj_u = 0;
548  double cell_l2error_adj_v = 0;
549  for(unsigned int iquad = 0; iquad < n_quad_pts; ++iquad){
550  std::fill(soln_at_q_u.begin(), soln_at_q_u.end(), 0);
551  std::fill(soln_at_q_v.begin(), soln_at_q_v.end(), 0);
552  std::fill(adj_at_q_u.begin(), adj_at_q_u.end(), 0);
553  std::fill(adj_at_q_v.begin(), adj_at_q_v.end(), 0);
554 
555  for (unsigned int idof = 0; idof < fe_values_extra.dofs_per_cell; ++idof) {
556  const unsigned int istate = fe_values_extra.get_fe().system_to_component_index(idof).first;
557  soln_at_q_u[istate] += dg_u->solution[dofs_indices[idof]] * fe_values_extra.shape_value_component(idof, iquad, istate);
558  soln_at_q_v[istate] += dg_v->solution[dofs_indices[idof]] * fe_values_extra.shape_value_component(idof, iquad, istate);
559  adj_at_q_u[istate] += adjoint_u[dofs_indices[idof]] * fe_values_extra.shape_value_component(idof, iquad, istate);
560  adj_at_q_v[istate] += adjoint_v[dofs_indices[idof]] * fe_values_extra.shape_value_component(idof, iquad, istate);
561  }
562 
563  const dealii::Point<dim,double> qpoint = (fe_values_extra.quadrature_point(iquad));
564 
565  for (int istate = 0; istate < nstate; ++istate){
566  const double soln_exact_u = physics_u_double->manufactured_solution_function->value(qpoint, istate);
567  const double soln_exact_v = physics_v_double->manufactured_solution_function->value(qpoint, istate);
568 
569  // comparing the converged solution to the manufactured solution
570  cell_l2error_soln_u += pow(soln_at_q_u[istate] - soln_exact_u, 2) * fe_values_extra.JxW(iquad);
571  cell_l2error_soln_v += pow(soln_at_q_v[istate] - soln_exact_v, 2) * fe_values_extra.JxW(iquad);
572 
573  // adjoint should convert to the manufactured solution of the opposing case
574  cell_l2error_adj_u += pow(adj_at_q_u[istate] - soln_exact_v, 2) * fe_values_extra.JxW(iquad);
575  cell_l2error_adj_v += pow(adj_at_q_v[istate] - soln_exact_u, 2) * fe_values_extra.JxW(iquad);
576 
577  // std::cout << "Adjoint value is = " << adj_at_q_u[istate] << std::endl << "and the exact value is = " << adj_exact_u << std::endl;
578  }
579  }
580 
581  // std::cout << "Cell value is (for u) = " << cell_l2error_adj_u << std::endl;
582 
583  // Storing the cell-wise error terms
584  cellError_soln_u[cell->active_cell_index()] = cell_l2error_soln_u;
585  cellError_soln_v[cell->active_cell_index()] = cell_l2error_soln_v;
586  cellError_adj_u[cell->active_cell_index()] = cell_l2error_adj_u;
587  cellError_adj_v[cell->active_cell_index()] = cell_l2error_adj_v;
588 
589  // Adding contributions to the global error
590  l2error_soln_u += cell_l2error_soln_u;
591  l2error_soln_v += cell_l2error_soln_v;
592  l2error_adj_u += cell_l2error_adj_u;
593  l2error_adj_v += cell_l2error_adj_v;
594  }
595  const double l2error_soln_u_mpi_sum = std::sqrt(dealii::Utilities::MPI::sum(l2error_soln_u, mpi_communicator));
596  const double l2error_soln_v_mpi_sum = std::sqrt(dealii::Utilities::MPI::sum(l2error_soln_v, mpi_communicator));
597  const double l2error_adj_u_mpi_sum = std::sqrt(dealii::Utilities::MPI::sum(l2error_adj_u, mpi_communicator));
598  const double l2error_adj_v_mpi_sum = std::sqrt(dealii::Utilities::MPI::sum(l2error_adj_v, mpi_communicator));
599 
600  // outputing the solutions obtained
601  if(dim == 1){
602  const std::string filename_u = "sol-u-" + std::to_string(igrid) + ".gnuplot";
603  std::ofstream gnuplot_output_u(filename_u);
604  dealii::DataOut<dim, dealii::DoFHandler<dim>> data_out_u;
605  data_out_u.attach_dof_handler(dg_u->dof_handler);
606  data_out_u.add_data_vector(dg_u->solution, "u");
607  data_out_u.build_patches();
608  data_out_u.write_gnuplot(gnuplot_output_u);
609 
610  const std::string filename_v = "sol-v-" + std::to_string(igrid) + ".gnuplot";
611  std::ofstream gnuplot_output_v(filename_v);
612  dealii::DataOut<dim, dealii::DoFHandler<dim>> data_out_v;
613  data_out_v.attach_dof_handler(dg_v->dof_handler);
614  data_out_v.add_data_vector(dg_v->solution, "u");
615  data_out_v.build_patches();
616  data_out_v.write_gnuplot(gnuplot_output_v);
617  }else{
618  pcout << "Outputting the results" << std::endl;
619  // // outputing the adjoint
620  // adj_u.output_results_vtk(10*poly_degree + igrid);
621  // adj_v.output_results_vtk(100 + 10*poly_degree + igrid);
622 
623  // // initializing vectors for source terms and manufactured solution
624  // dealii::LinearAlgebra::distributed::Vector<double> source_u; source_u.reinit(dg_u->solution);
625  // dealii::LinearAlgebra::distributed::Vector<double> source_v; source_v.reinit(dg_u->solution);
626  // dealii::LinearAlgebra::distributed::Vector<double> manu_u; manu_u.reinit(dg_u->solution);
627  // dealii::LinearAlgebra::distributed::Vector<double> manu_v; manu_v.reinit(dg_u->solution);
628  // need to add a dof-wise loop to output these quantities. Haven't figured out how yet that doesn't rely on interpolation of a dealii::Function
629 
630  // setting up dataout and outputing results
631  dealii::DataOut<dim, dealii::DoFHandler<dim>> data_out;
632  data_out.attach_dof_handler(dg_u->dof_handler);
633 
634  // // can't use this post processor as it gives them the same name
635  // const std::unique_ptr< dealii::DataPostprocessor<dim> > post_processor_u = Postprocess::PostprocessorFactory<dim>::create_Postprocessor(all_parameters);
636  // const std::unique_ptr< dealii::DataPostprocessor<dim> > post_processor_v = Postprocess::PostprocessorFactory<dim>::create_Postprocessor(all_parameters);
637  // data_out.add_data_vector(dg_u->solution, *post_processor_u);
638  // data_out.add_data_vector(dg_v->solution, *post_processor_v);
639 
640  dealii::Vector<float> subdomain(dg_u->triangulation->n_active_cells());
641  for (unsigned int i = 0; i < subdomain.size(); ++i) {
642  subdomain(i) = dg_u->triangulation->locally_owned_subdomain();
643  }
644  data_out.add_data_vector(subdomain, "subdomain", dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
645 
646  std::vector<unsigned int> active_fe_indices;
647  dg_u->dof_handler.get_active_fe_indices(active_fe_indices);
648  dealii::Vector<double> active_fe_indices_dealiivector(active_fe_indices.begin(), active_fe_indices.end());
649  dealii::Vector<double> cell_poly_degree = active_fe_indices_dealiivector;
650 
651  data_out.add_data_vector(active_fe_indices_dealiivector, "PolynomialDegree", dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
652 
653  std::vector<std::string> solution_names_u;
654  std::vector<std::string> solution_names_v;
655  std::vector<std::string> residual_names_u;
656  std::vector<std::string> residual_names_v;
657  std::vector<std::string> dIdw_names_u;
658  std::vector<std::string> dIdw_names_v;
659  std::vector<std::string> adjoint_names_u;
660  std::vector<std::string> adjoint_names_v;
661  for(int s=0;s<nstate;++s) {
662  std::string varname0_u = "state" + dealii::Utilities::int_to_string(s,1) + "_u";
663  std::string varname0_v = "state" + dealii::Utilities::int_to_string(s,1) + "_v";
664  std::string varname1_u = "residual" + dealii::Utilities::int_to_string(s,1) + "_u";
665  std::string varname1_v = "residual" + dealii::Utilities::int_to_string(s,1) + "_v";
666  std::string varname2_u = "dIdw" + dealii::Utilities::int_to_string(s,1) + "_u";
667  std::string varname2_v = "dIdw" + dealii::Utilities::int_to_string(s,1) + "_v";
668  std::string varname3_u = "psi" + dealii::Utilities::int_to_string(s,1) + "_u";
669  std::string varname3_v = "psi" + dealii::Utilities::int_to_string(s,1) + "_v";
670 
671  solution_names_u.push_back(varname0_u);
672  solution_names_v.push_back(varname0_v);
673  residual_names_u.push_back(varname1_u);
674  residual_names_v.push_back(varname1_v);
675  dIdw_names_u.push_back(varname2_u);
676  dIdw_names_v.push_back(varname2_v);
677  adjoint_names_u.push_back(varname3_u);
678  adjoint_names_v.push_back(varname3_v);
679  }
680 
681  data_out.add_data_vector(dg_u->solution, solution_names_u, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
682  data_out.add_data_vector(dg_v->solution, solution_names_v, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
683  data_out.add_data_vector(dg_u->right_hand_side, residual_names_u, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
684  data_out.add_data_vector(dg_v->right_hand_side, residual_names_v, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
685  data_out.add_data_vector(adj_u.dIdw_coarse, dIdw_names_u, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
686  data_out.add_data_vector(adj_v.dIdw_coarse, dIdw_names_v, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
687  data_out.add_data_vector(adj_u.adjoint_coarse, adjoint_names_u, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
688  data_out.add_data_vector(adj_v.adjoint_coarse, adjoint_names_v, dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_dof_data);
689 
690  data_out.add_data_vector(cellError_soln_u, "soln_u_err", dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
691  data_out.add_data_vector(cellError_soln_v, "soln_v_err", dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
692  data_out.add_data_vector(cellError_adj_u, "adj_u_err", dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
693  data_out.add_data_vector(cellError_adj_v, "adj_v_err", dealii::DataOut_DoFData<dealii::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
694 
695  const int iproc = dealii::Utilities::MPI::this_mpi_process(mpi_communicator);
696  data_out.build_patches();
697 
698  // creating the files
699  std::string filename = "diffusion_exact_adjoint-" ;
700  filename += dealii::Utilities::int_to_string(dim, 1) + "D-";
701  filename += dealii::Utilities::int_to_string(poly_degree, 1) + "P-";
702  filename += dealii::Utilities::int_to_string(igrid, 4) + ".";
703  filename += dealii::Utilities::int_to_string(iproc, 4);
704  filename += ".vtu";
705  std::ofstream output(filename);
706  data_out.write_vtu(output);
707 
708  if (iproc == 0) {
709  std::vector<std::string> filenames;
710  for (unsigned int iproc = 0; iproc < dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); ++iproc) {
711  std::string fn = "diffusion_exact_adjoint-";
712  fn += dealii::Utilities::int_to_string(dim, 1) + "D-";
713  fn += dealii::Utilities::int_to_string(poly_degree, 1) + "P-";
714  fn += dealii::Utilities::int_to_string(igrid, 4) + ".";
715  fn += dealii::Utilities::int_to_string(iproc, 4);
716  fn += ".vtu";
717  filenames.push_back(fn);
718  }
719  std::string master_fn = "diffusion_exact_adjoint-";
720  master_fn += dealii::Utilities::int_to_string(dim, 1) +"D-";
721  master_fn += dealii::Utilities::int_to_string(poly_degree, 1) +"P-";
722  master_fn += dealii::Utilities::int_to_string(igrid, 4) + ".pvtu";
723  std::ofstream master_output(master_fn);
724  data_out.write_pvtu_record(master_output, filenames);
725  }
726  }
727 
728  // adding terms to the convergence table
729  const int n_dofs = dg_u->dof_handler.n_dofs();
730  const double dx = 1.0/pow(n_dofs,(1.0/dim));
731 
732  // adding terms to the table
733  convergence_table.add_value("p", poly_degree);
734  convergence_table.add_value("cells", grid->n_global_active_cells());
735  convergence_table.add_value("DoFs", n_dofs);
736  convergence_table.add_value("dx", dx);
737  convergence_table.add_value("soln_u_val", functional_val_u);
738  convergence_table.add_value("soln_v_val", functional_val_v);
739  convergence_table.add_value("soln_u_err", error_functional_u);
740  convergence_table.add_value("soln_v_err", error_functional_v);
741  convergence_table.add_value("soln_u_L2_err", l2error_soln_u_mpi_sum);
742  convergence_table.add_value("soln_v_L2_err", l2error_soln_v_mpi_sum);
743  convergence_table.add_value("adj_u_L2_err", l2error_adj_u_mpi_sum);
744  convergence_table.add_value("adj_v_L2_err", l2error_adj_v_mpi_sum);
745 
746  // storing in vectors for convergence checing
747  grid_size[igrid] = dx;
748  output_error_u[igrid] = error_functional_u;
749  output_error_v[igrid] = error_functional_v;
750  soln_error_u[igrid] = l2error_soln_u_mpi_sum;
751  soln_error_v[igrid] = l2error_soln_v_mpi_sum;
752  adj_error_u[igrid] = l2error_adj_u_mpi_sum;
753  adj_error_v[igrid] = l2error_adj_v_mpi_sum;
754  }
755 
756  // obtaining the convergence rates
757  convergence_table.evaluate_convergence_rates("soln_u_err", "cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
758  convergence_table.evaluate_convergence_rates("soln_v_err", "cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
759  convergence_table.evaluate_convergence_rates("soln_u_L2_err", "cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
760  convergence_table.evaluate_convergence_rates("soln_v_L2_err", "cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
761  convergence_table.evaluate_convergence_rates("adj_u_L2_err", "cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
762  convergence_table.evaluate_convergence_rates("adj_v_L2_err", "cells", dealii::ConvergenceTable::reduction_rate_log2, dim);
763  convergence_table.set_scientific("dx",true);
764  convergence_table.set_scientific("soln_u_val",true);
765  convergence_table.set_scientific("soln_v_val",true);
766  convergence_table.set_scientific("soln_u_err",true);
767  convergence_table.set_scientific("soln_v_err",true);
768  convergence_table.set_scientific("soln_u_L2_err",true);
769  convergence_table.set_scientific("soln_v_L2_err",true);
770  convergence_table.set_scientific("adj_u_L2_err",true);
771  convergence_table.set_scientific("adj_v_L2_err",true);
772 
773  // adding it to the final list
774  convergence_table_vector.push_back(convergence_table);
775 
776  // setting slope targets for convergence orders
777  // const double expected_slope_error_functional = 2*poly_degree + 1;
778  // const double expected_slope_l2error_soln = poly_degree + 1;
779  const double expected_slope_l2error_adj = poly_degree + 1;
780 
781  // evaluating the average slopes from the last two steps
782  // const double avg_slope_error_functional_u = eval_avg_slope(output_error_u, grid_size, n_grids);
783  // const double avg_slope_error_functional_v = eval_avg_slope(output_error_v, grid_size, n_grids);
784  // const double avg_slope_error_l2error_soln_u = eval_avg_slope( soln_error_u, grid_size, n_grids);
785  // const double avg_slope_error_l2error_soln_v = eval_avg_slope( soln_error_v, grid_size, n_grids);
786  const double avg_slope_error_l2error_adj_u = eval_avg_slope( adj_error_u, grid_size, n_grids);
787  const double avg_slope_error_l2error_adj_v = eval_avg_slope( adj_error_v, grid_size, n_grids);
788 
789  // diffference from the expected
790  // const double diff_slope_error_functional_u = avg_slope_error_functional_u - expected_slope_error_functional;
791  // const double diff_slope_error_functional_v = avg_slope_error_functional_v - expected_slope_error_functional;
792  // const double diff_slope_error_l2error_soln_u = avg_slope_error_l2error_soln_u - expected_slope_l2error_soln;
793  // const double diff_slope_error_l2error_soln_v = avg_slope_error_l2error_soln_v - expected_slope_l2error_soln;
794  const double diff_slope_error_l2error_adj_u = avg_slope_error_l2error_adj_u - expected_slope_l2error_adj;
795  const double diff_slope_error_l2error_adj_v = avg_slope_error_l2error_adj_v - expected_slope_l2error_adj;
796 
797  // tolerance set from the input file
798  double slope_deficit_tolerance = -std::abs(manu_grid_conv_param.slope_deficit_tolerance);
799 
800  // // performing the actual checks
801  // if(diff_slope_error_functional_u < slope_deficit_tolerance){
802  // pcout << "Convergence order not achieved for functional_u." << std::endl
803  // << "Average order of " << avg_slope_error_functional_u << " instead of expected " << expected_slope_error_functional << std::endl;
804  // n_fail++;
805  // }
806  // if(diff_slope_error_functional_v < slope_deficit_tolerance){
807  // pcout << "Convergence order not achieved for functional_v." << std::endl
808  // << "Average order of " << avg_slope_error_functional_v << " instead of expected " << expected_slope_error_functional << std::endl;
809  // n_fail++;
810  // }
811  // if(diff_slope_error_l2error_soln_u < slope_deficit_tolerance){
812  // pcout << "Convergence order not achieved for l2error_soln_u." << std::endl
813  // << "Average order of " << avg_slope_error_l2error_soln_u << " instead of expected " << expected_slope_l2error_soln << std::endl;
814  // n_fail++;
815  // }
816  // if(diff_slope_error_l2error_soln_v < slope_deficit_tolerance){
817  // pcout << "Convergence order not achieved for l2error_soln_v." << std::endl
818  // << "Average order of " << avg_slope_error_l2error_soln_v << " instead of expected " << expected_slope_l2error_soln << std::endl;
819  // n_fail++;
820  // }
821  if(diff_slope_error_l2error_adj_u < slope_deficit_tolerance){
822  pcout << "Convergence order not achieved for l2error_adj_u." << std::endl
823  << "Average order of " << avg_slope_error_l2error_adj_u << " instead of expected " << expected_slope_l2error_adj << std::endl;
824  n_fail++;
825  }
826  if(diff_slope_error_l2error_adj_v < slope_deficit_tolerance){
827  pcout << "Convergence order not achieved for l2error_adj_v." << std::endl
828  << "Average order of " << avg_slope_error_l2error_adj_v << " instead of expected " << expected_slope_l2error_adj << std::endl;
829  n_fail++;
830  }
831  }
832 
833  pcout << std::endl << std::endl << std::endl << std::endl;
834  pcout << " ********************************************" << std::endl;
835  pcout << " Convergence summary" << std::endl;
836  pcout << " ********************************************" << std::endl;
837  for (auto conv = convergence_table_vector.begin(); conv!=convergence_table_vector.end(); conv++) {
838  if (pcout.is_active()) conv->write_text(pcout.get_stream());
839  pcout << " ********************************************" << std::endl;
840  }
841 
842  return n_fail;
843 }
844 
845 // computes the average error of the last two slopes
846 double eval_avg_slope(std::vector<double> error, std::vector<double> grid_size, unsigned int n_grids){
847  const double last_slope_error = log(error[n_grids-1]/error[n_grids-2])/(log(grid_size[n_grids-1]/grid_size[n_grids-2]));
848  double prev_slope_error = last_slope_error;
849  if(n_grids > 2){
850  prev_slope_error = log(error[n_grids-2]/error[n_grids-3])/(log(grid_size[n_grids-2]/grid_size[n_grids-3]));
851  }
852  return 0.5*(last_slope_error+prev_slope_error);
853 }
854 
856 
857 } // Tests namespace
858 } // PHiLiP namespace
real2 evaluate_volume_integrand(const PHiLiP::Physics::PhysicsBase< dim, nstate, real2 > &physics, const dealii::Point< dim, real2 > &phys_coord, const std::array< real2, nstate > &soln_at_q, const std::array< dealii::Tensor< 1, dim, real2 >, nstate > &soln_grad_at_q) const
Templated volume integrand.
real objective_function(const dealii::Point< dim, real > &pos) const override
objective function = f
unsigned int degree_start
First polynomial degree to start the loop. If diffusion, must be at least 1.
dealii::LinearAlgebra::distributed::Vector< real > dIdw_coarse
functional derivative (on the coarse grid)
Definition: adjoint.h:134
PartialDifferentialEquation pde_type
Store the PDE type to be solved.
dealii::LinearAlgebra::distributed::Vector< real > adjoint_coarse
coarse grid adjoint ( )
Definition: adjoint.h:138
virtual real objective_function(const dealii::Point< dim, real > &pos) const =0
defnined directly as part of the physics to make passing to the functional simpler ...
Parameters related to the manufactured convergence study.
Base class from which Advection, Diffusion, ConvectionDiffusion, and Euler is derived.
Definition: physics.h:34
Adjoint class.
Definition: adjoint.h:39
const MPI_Comm mpi_communicator
MPI communicator.
Definition: tests.h:39
PartialDifferentialEquation
Possible Partial Differential Equations to solve.
std::array< real, nstate > source_term(const dealii::Point< dim, real > &pos, const std::array< real, nstate > &solution, const real current_time, const dealii::types::global_dof_index cell_index) const override
source term = f
Files for the baseline physics.
Definition: ADTypes.hpp:10
DiffusionExactAdjoint(const Parameters::AllParameters *const parameters_input)
Constructor to call the TestsBase constructor to set parameters = parameters_input.
ManufacturedSolutionParam manufactured_solution_param
Associated manufactured solution parameters.
unsigned int dimension
Number of dimensions. Note that it has to match the executable PHiLiP_xD.
Main parameter class that contains the various other sub-parameter classes.
real objective_function(const dealii::Point< dim, real > &pos) const override
objective function = g
ManufacturedConvergenceStudyParam manufactured_convergence_study_param
Contains parameters for manufactured convergence study.
dealii::LinearAlgebra::distributed::Vector< real > coarse_grid_adjoint()
Computes the coarse grid adjoint.
Definition: adjoint.cpp:187
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &pos, const unsigned int istate=0) const override
Gradient of the manufactured solution.
const Parameters::AllParameters *const all_parameters
Pointer to all parameters.
Definition: tests.h:20
static std::shared_ptr< DGBase< dim, real, MeshType > > create_discontinuous_galerkin(const Parameters::AllParameters *const parameters_input, const unsigned int degree, const unsigned int max_degree_input, const unsigned int grid_degree_input, const std::shared_ptr< Triangulation > triangulation_input)
Creates a derived object DG, but returns it as DGBase.
Definition: dg_factory.cpp:10
real value(const dealii::Point< dim, real > &pos, const unsigned int istate=0) const override
overriding the function for the value and gradient
std::array< real, nstate > source_term(const dealii::Point< dim, real > &pos, const std::array< real, nstate > &solution, const real current_time, const dealii::types::global_dof_index cell_index) const override
source term = g
bool use_manufactured_source_term
Uses non-zero source term based on the manufactured solution and the PDE.
Abstract class templated on the number of state variables.
dealii::ConditionalOStream pcout
ConditionalOStream.
Definition: tests.h:45
dealii::Tensor< 1, dim, real > gradient(const dealii::Point< dim, real > &pos, const unsigned int istate=0) const override
Gradient of the manufactured solution.
std::vector< int > get_number_1d_cells(const int ngrids) const
Evaluates the number of cells to generate the grids for 1D grid based on input file.
Definition: tests.cpp:71
real value(const dealii::Point< dim, real > &pos, const unsigned int istate=0) const override
overriding the function for the value and gradient
static std::shared_ptr< ODESolverBase< dim, real, MeshType > > create_ODESolver(std::shared_ptr< DGBase< dim, real, MeshType > > dg_input)
Creates either implicit or explicit ODE solver based on parameter value(no POD basis given) ...
void set_physics(std::shared_ptr< Physics::PhysicsBase< dim, nstate, real > > pde_physics_double_input, std::shared_ptr< Physics::PhysicsBase< dim, nstate, FadType > > pde_physics_fad_input, std::shared_ptr< Physics::PhysicsBase< dim, nstate, RadType > > pde_physics_rad_input, std::shared_ptr< Physics::PhysicsBase< dim, nstate, FadFadType > > pde_physics_fad_fad_input, std::shared_ptr< Physics::PhysicsBase< dim, nstate, RadFadType > > pde_physics_rad_fad_input)
Base class of all the tests.
Definition: tests.h:17
parent class to add the objective function directly to physics as a virtual class ...