[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
physics.cpp
1 #include <assert.h>
2 #include <cmath>
3 #include <vector>
4 #include <stdlib.h>
5 #include <deal.II/base/utilities.h>
6 #include <deal.II/base/mpi.h>
7 
8 #include "ADTypes.hpp"
9 
10 #include "physics.h"
11 
12 namespace PHiLiP {
13 namespace Physics {
14 
15 template <int dim, int nstate, typename real>
16 PhysicsBase<dim,nstate,real>::PhysicsBase(
17  const Parameters::AllParameters *const parameters_input,
18  const bool has_nonzero_diffusion_input,
19  const bool has_nonzero_physical_source_input,
20  const dealii::Tensor<2,3,double> input_diffusion_tensor,
21  std::shared_ptr< ManufacturedSolutionFunction<dim,real> > manufactured_solution_function_input)
22  : has_nonzero_diffusion(has_nonzero_diffusion_input)
23  , has_nonzero_physical_source(has_nonzero_physical_source_input)
24  , all_parameters(parameters_input)
25  , non_physical_behavior_type(all_parameters->non_physical_behavior_type)
26  , manufactured_solution_function(manufactured_solution_function_input)
27  , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0)
28 {
29  // if provided with a null ptr, give it the default manufactured solution
30  // currently only necessary for the unit test
31  if(!manufactured_solution_function)
32  manufactured_solution_function = std::make_shared<ManufacturedSolutionSine<dim,real>>(nstate);
33 
34  // anisotropic diffusion matrix
35  diffusion_tensor[0][0] = input_diffusion_tensor[0][0];
36  if constexpr(dim >= 2) {
37  diffusion_tensor[0][1] = input_diffusion_tensor[0][1];
38  diffusion_tensor[1][0] = input_diffusion_tensor[1][0];
39  diffusion_tensor[1][1] = input_diffusion_tensor[1][1];
40  }
41  if constexpr(dim >= 3) {
42  diffusion_tensor[0][2] = input_diffusion_tensor[0][2];
43  diffusion_tensor[2][0] = input_diffusion_tensor[2][0];
44  diffusion_tensor[1][2] = input_diffusion_tensor[1][2];
45  diffusion_tensor[2][1] = input_diffusion_tensor[2][1];
46  diffusion_tensor[2][2] = input_diffusion_tensor[2][2];
47  }
48 }
49 
50 template <int dim, int nstate, typename real>
51 PhysicsBase<dim,nstate,real>::PhysicsBase(
52  const Parameters::AllParameters *const parameters_input,
53  const bool has_nonzero_diffusion_input,
54  const bool has_nonzero_physical_source_input,
55  std::shared_ptr< ManufacturedSolutionFunction<dim,real> > manufactured_solution_function_input)
56  : PhysicsBase<dim,nstate,real>(
57  parameters_input,
58  has_nonzero_diffusion_input,
59  has_nonzero_physical_source_input,
60  Parameters::ManufacturedSolutionParam::get_default_diffusion_tensor(),
61  manufactured_solution_function_input)
62 { }
63 
64 template <int dim, int nstate, typename real>
65 std::array<dealii::Tensor<1,dim,real>,nstate> PhysicsBase<dim,nstate,real>::convective_numerical_split_flux (
66  const std::array<real,nstate> &/*conservative_soln1*/,
67  const std::array<real,nstate> &/*conservative_soln2*/) const
68 {
69  pcout << "ERROR: convective_numerical_split_flux() has not yet been implemented for (overridden by) the selected PDE. Aborting..." <<std::flush;
70  std::abort();
71  std::array<dealii::Tensor<1,dim,real>,nstate> dummy;
72  return dummy;
73 }
74 
75 template <int dim, int nstate, typename real>
76 real PhysicsBase<dim,nstate,real>
77 ::max_convective_normal_eigenvalue (
78  const std::array<real,nstate> &conservative_soln,
79  const dealii::Tensor<1,dim,real> &/*normal*/) const
80 {
81  return max_convective_eigenvalue(conservative_soln);
82 }
83 
84 /*
85 template <int dim, int nstate, typename real>
86 std::array<dealii::Tensor<1,dim,real>,nstate> PhysicsBase<dim,nstate,real>
87 ::artificial_dissipative_flux (
88  const real viscosity_coefficient,
89  const std::array<real,nstate> &,//solution,
90  const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient) const
91 {
92  std::array<dealii::Tensor<1,dim,real>,nstate> diss_flux;
93  for (int i=0; i<nstate; i++) {
94  for (int d=0; d<dim; d++) {
95  diss_flux[i][d] = -viscosity_coefficient*(solution_gradient[i][d]);
96  }
97  }
98  return diss_flux;
99 }
100 */
101 template <int dim, int nstate, typename real>
102 std::array<real,nstate> PhysicsBase<dim,nstate,real>
103 ::artificial_source_term (
104  const real viscosity_coefficient,
105  const dealii::Point<dim,real> &pos,
106  const std::array<real,nstate> &/*solution*/) const
107 {
108  std::array<real,nstate> source;
109 
110  dealii::Tensor<2,dim,double> artificial_diffusion_tensor;
111  for (int i=0;i<dim;i++)
112  for (int j=0;j<dim;j++)
113  artificial_diffusion_tensor[i][j] = (i==j) ? 1.0 : 0.0;
114 
115  for (int istate=0; istate<nstate; istate++) {
116  dealii::SymmetricTensor<2,dim,real> manufactured_hessian = this->manufactured_solution_function->hessian (pos, istate);
117  //source[istate] = -viscosity_coefficient*scalar_product(artificial_diffusion_tensor,manufactured_hessian);
118  source[istate] = 0.0;
119  for (int dr=0; dr<dim; ++dr) {
120  for (int dc=0; dc<dim; ++dc) {
121  source[istate] += artificial_diffusion_tensor[dr][dc] * manufactured_hessian[dr][dc];
122  }
123  }
124  source[istate] *= -viscosity_coefficient;
125  }
126  return source;
127 }
128 
129 template <int dim, int nstate, typename real>
130 std::array<real,nstate> PhysicsBase<dim,nstate,real>
131 ::physical_source_term (
132  const dealii::Point<dim,real> &/*pos*/,
133  const std::array<real,nstate> &/*solution*/,
134  const std::array<dealii::Tensor<1,dim,real>,nstate> &/*solution_gradient*/,
135  const dealii::types::global_dof_index /*cell_index*/) const
136 {
137  std::array<real,nstate> physical_source;
138  for (int i=0; i<nstate; i++) {
139  physical_source[i] = 0;
140  }
141  return physical_source;
142 }
143 
144 template <int dim, int nstate, typename real>
145 dealii::Vector<double> PhysicsBase<dim,nstate,real>::post_compute_derived_quantities_vector (
146  const dealii::Vector<double> &uh,
147  const std::vector<dealii::Tensor<1,dim> > &/*duh*/,
148  const std::vector<dealii::Tensor<2,dim> > &/*dduh*/,
149  const dealii::Tensor<1,dim> &/*normals*/,
150  const dealii::Point<dim> &/*evaluation_points*/) const
151 {
152  dealii::Vector<double> computed_quantities(nstate);
153  for (unsigned int s=0; s<nstate; ++s) {
154  computed_quantities(s) = uh(s);
155  }
156  return computed_quantities;
157 }
158 
159 template <int dim, int nstate, typename real>
160 dealii::Vector<double> PhysicsBase<dim,nstate,real>::post_compute_derived_quantities_scalar (
161  const double &uh,
162  const dealii::Tensor<1,dim> &/*duh*/,
163  const dealii::Tensor<2,dim> &/*dduh*/,
164  const dealii::Tensor<1,dim> &/*normals*/,
165  const dealii::Point<dim> &/*evaluation_points*/) const
166 {
167  assert(nstate == 1);
168  dealii::Vector<double> computed_quantities(nstate);
169  for (unsigned int s=0; s<nstate; ++s) {
170  computed_quantities(s) = uh;
171  }
172  return computed_quantities;
173 }
174 
175 template <int dim, int nstate, typename real>
176 std::vector<std::string> PhysicsBase<dim,nstate,real> ::post_get_names () const
177 {
178  std::vector<std::string> names;
179  for (unsigned int s=0; s<nstate; ++s) {
180  std::string varname = "state" + dealii::Utilities::int_to_string(s,1);
181  names.push_back(varname);
182  }
183  return names;
184 }
185 
186 template <int dim, int nstate, typename real>
187 std::vector<dealii::DataComponentInterpretation::DataComponentInterpretation> PhysicsBase<dim,nstate,real>
188 ::post_get_data_component_interpretation () const
189 {
190  namespace DCI = dealii::DataComponentInterpretation;
191  std::vector<DCI::DataComponentInterpretation> interpretation;
192  for (unsigned int s=0; s<nstate; ++s) {
193  interpretation.push_back (DCI::component_is_scalar);
194  }
195  return interpretation;
196 }
197 
198 template <int dim, int nstate, typename real>
199 dealii::UpdateFlags PhysicsBase<dim,nstate,real>
200 ::post_get_needed_update_flags () const
201 {
202  return dealii::update_values;
203 }
204 
205 template <int dim, int nstate, typename real>
206 template<typename real2>
207 real2 PhysicsBase<dim,nstate,real>
208 ::handle_non_physical_result(const std::string message) const
209 {
210  if (this->non_physical_behavior_type == NonPhysicalBehaviorEnum::abort_run) {
211  std::cout << "ERROR: Non-physical result has been detected. ";
212  if (!message.empty()) {
213  std::cout << std::endl << " Message: " << message << std::endl;
214  }
215  std::cout << " Aborting... " << std::endl << std::flush;
216  std::abort();
217  } else if (this->non_physical_behavior_type == NonPhysicalBehaviorEnum::print_warning) {
218  std::cout << "WARNING: Non-physical result has been detected at a node." << std::endl;
219  if (!message.empty()) {
220  std::cout << std::endl << " Message: " << message << std::endl;
221  }
222  } else if (this->non_physical_behavior_type == NonPhysicalBehaviorEnum::return_big_number) {
223  // do nothing -- assume that the test or iterative solver can handle this.
224  }
225 
226  return (real2)BIG_NUMBER;
227 }
228 
229 template class PhysicsBase < PHILIP_DIM, 1, double >;
230 template class PhysicsBase < PHILIP_DIM, 2, double >;
231 template class PhysicsBase < PHILIP_DIM, 3, double >;
232 template class PhysicsBase < PHILIP_DIM, 4, double >;
233 template class PhysicsBase < PHILIP_DIM, 5, double >;
234 template class PhysicsBase < PHILIP_DIM, 6, double >;
235 template class PhysicsBase < PHILIP_DIM, 8, double >;
236 
237 template class PhysicsBase < PHILIP_DIM, 1, FadType >;
238 template class PhysicsBase < PHILIP_DIM, 2, FadType >;
239 template class PhysicsBase < PHILIP_DIM, 3, FadType >;
240 template class PhysicsBase < PHILIP_DIM, 4, FadType >;
241 template class PhysicsBase < PHILIP_DIM, 5, FadType >;
242 template class PhysicsBase < PHILIP_DIM, 6, FadType >;
243 template class PhysicsBase < PHILIP_DIM, 8, FadType >;
244 
245 template class PhysicsBase < PHILIP_DIM, 1, RadType >;
246 template class PhysicsBase < PHILIP_DIM, 2, RadType >;
247 template class PhysicsBase < PHILIP_DIM, 3, RadType >;
248 template class PhysicsBase < PHILIP_DIM, 4, RadType >;
249 template class PhysicsBase < PHILIP_DIM, 5, RadType >;
250 template class PhysicsBase < PHILIP_DIM, 6, RadType >;
251 template class PhysicsBase < PHILIP_DIM, 8, RadType >;
252 
253 template class PhysicsBase < PHILIP_DIM, 1, FadFadType >;
254 template class PhysicsBase < PHILIP_DIM, 2, FadFadType >;
255 template class PhysicsBase < PHILIP_DIM, 3, FadFadType >;
256 template class PhysicsBase < PHILIP_DIM, 4, FadFadType >;
257 template class PhysicsBase < PHILIP_DIM, 5, FadFadType >;
258 template class PhysicsBase < PHILIP_DIM, 6, FadFadType >;
259 template class PhysicsBase < PHILIP_DIM, 8, FadFadType >;
260 
261 template class PhysicsBase < PHILIP_DIM, 1, RadFadType >;
262 template class PhysicsBase < PHILIP_DIM, 2, RadFadType >;
263 template class PhysicsBase < PHILIP_DIM, 3, RadFadType >;
264 template class PhysicsBase < PHILIP_DIM, 4, RadFadType >;
265 template class PhysicsBase < PHILIP_DIM, 5, RadFadType >;
266 template class PhysicsBase < PHILIP_DIM, 6, RadFadType >;
267 template class PhysicsBase < PHILIP_DIM, 8, RadFadType >;
268 
269 //==============================================================================
270 // -> Templated member functions: // could be automated later on using Boost MPL
271 //------------------------------------------------------------------------------
272 // -- handle_non_physical_result
273 template double PhysicsBase < PHILIP_DIM, 1, double >::handle_non_physical_result<double>(const std::string message) const;
274 template double PhysicsBase < PHILIP_DIM, 2, double >::handle_non_physical_result<double>(const std::string message) const;
275 template double PhysicsBase < PHILIP_DIM, 3, double >::handle_non_physical_result<double>(const std::string message) const;
276 template double PhysicsBase < PHILIP_DIM, 4, double >::handle_non_physical_result<double>(const std::string message) const;
277 template double PhysicsBase < PHILIP_DIM, 5, double >::handle_non_physical_result<double>(const std::string message) const;
278 template double PhysicsBase < PHILIP_DIM, 6, double >::handle_non_physical_result<double>(const std::string message) const;
279 template double PhysicsBase < PHILIP_DIM, 8, double >::handle_non_physical_result<double>(const std::string message) const;
280 
281 template FadType PhysicsBase < PHILIP_DIM, 1, FadType >::handle_non_physical_result<FadType>(const std::string message) const;
282 template FadType PhysicsBase < PHILIP_DIM, 2, FadType >::handle_non_physical_result<FadType>(const std::string message) const;
283 template FadType PhysicsBase < PHILIP_DIM, 3, FadType >::handle_non_physical_result<FadType>(const std::string message) const;
284 template FadType PhysicsBase < PHILIP_DIM, 4, FadType >::handle_non_physical_result<FadType>(const std::string message) const;
285 template FadType PhysicsBase < PHILIP_DIM, 5, FadType >::handle_non_physical_result<FadType>(const std::string message) const;
286 template FadType PhysicsBase < PHILIP_DIM, 6, FadType >::handle_non_physical_result<FadType>(const std::string message) const;
287 template FadType PhysicsBase < PHILIP_DIM, 8, FadType >::handle_non_physical_result<FadType>(const std::string message) const;
288 
289 template RadType PhysicsBase < PHILIP_DIM, 1, RadType >::handle_non_physical_result<RadType>(const std::string message) const;
290 template RadType PhysicsBase < PHILIP_DIM, 2, RadType >::handle_non_physical_result<RadType>(const std::string message) const;
291 template RadType PhysicsBase < PHILIP_DIM, 3, RadType >::handle_non_physical_result<RadType>(const std::string message) const;
292 template RadType PhysicsBase < PHILIP_DIM, 4, RadType >::handle_non_physical_result<RadType>(const std::string message) const;
293 template RadType PhysicsBase < PHILIP_DIM, 5, RadType >::handle_non_physical_result<RadType>(const std::string message) const;
294 template RadType PhysicsBase < PHILIP_DIM, 6, RadType >::handle_non_physical_result<RadType>(const std::string message) const;
295 template RadType PhysicsBase < PHILIP_DIM, 8, RadType >::handle_non_physical_result<RadType>(const std::string message) const;
296 
297 template FadFadType PhysicsBase < PHILIP_DIM, 1, FadFadType >::handle_non_physical_result<FadFadType>(const std::string message) const;
298 template FadFadType PhysicsBase < PHILIP_DIM, 2, FadFadType >::handle_non_physical_result<FadFadType>(const std::string message) const;
299 template FadFadType PhysicsBase < PHILIP_DIM, 3, FadFadType >::handle_non_physical_result<FadFadType>(const std::string message) const;
300 template FadFadType PhysicsBase < PHILIP_DIM, 4, FadFadType >::handle_non_physical_result<FadFadType>(const std::string message) const;
301 template FadFadType PhysicsBase < PHILIP_DIM, 5, FadFadType >::handle_non_physical_result<FadFadType>(const std::string message) const;
302 template FadFadType PhysicsBase < PHILIP_DIM, 6, FadFadType >::handle_non_physical_result<FadFadType>(const std::string message) const;
303 template FadFadType PhysicsBase < PHILIP_DIM, 8, FadFadType >::handle_non_physical_result<FadFadType>(const std::string message) const;
304 
305 template RadFadType PhysicsBase < PHILIP_DIM, 1, RadFadType >::handle_non_physical_result<RadFadType>(const std::string message) const;
306 template RadFadType PhysicsBase < PHILIP_DIM, 2, RadFadType >::handle_non_physical_result<RadFadType>(const std::string message) const;
307 template RadFadType PhysicsBase < PHILIP_DIM, 3, RadFadType >::handle_non_physical_result<RadFadType>(const std::string message) const;
308 template RadFadType PhysicsBase < PHILIP_DIM, 4, RadFadType >::handle_non_physical_result<RadFadType>(const std::string message) const;
309 template RadFadType PhysicsBase < PHILIP_DIM, 5, RadFadType >::handle_non_physical_result<RadFadType>(const std::string message) const;
310 template RadFadType PhysicsBase < PHILIP_DIM, 6, RadFadType >::handle_non_physical_result<RadFadType>(const std::string message) const;
311 template RadFadType PhysicsBase < PHILIP_DIM, 8, RadFadType >::handle_non_physical_result<RadFadType>(const std::string message) const;
312  // -- -- instantiate all the real types with real2 = FadType for automatic differentiation in NavierStokes::dissipative_flux_directional_jacobian()
313 template FadType PhysicsBase < PHILIP_DIM, 1, double >::handle_non_physical_result<FadType>(const std::string message) const;
314 template FadType PhysicsBase < PHILIP_DIM, 2, double >::handle_non_physical_result<FadType>(const std::string message) const;
315 template FadType PhysicsBase < PHILIP_DIM, 3, double >::handle_non_physical_result<FadType>(const std::string message) const;
316 template FadType PhysicsBase < PHILIP_DIM, 4, double >::handle_non_physical_result<FadType>(const std::string message) const;
317 template FadType PhysicsBase < PHILIP_DIM, 5, double >::handle_non_physical_result<FadType>(const std::string message) const;
318 template FadType PhysicsBase < PHILIP_DIM, 6, double >::handle_non_physical_result<FadType>(const std::string message) const;
319 template FadType PhysicsBase < PHILIP_DIM, 8, double >::handle_non_physical_result<FadType>(const std::string message) const;
320 
321 template FadType PhysicsBase < PHILIP_DIM, 1, RadType >::handle_non_physical_result<FadType>(const std::string message) const;
322 template FadType PhysicsBase < PHILIP_DIM, 2, RadType >::handle_non_physical_result<FadType>(const std::string message) const;
323 template FadType PhysicsBase < PHILIP_DIM, 3, RadType >::handle_non_physical_result<FadType>(const std::string message) const;
324 template FadType PhysicsBase < PHILIP_DIM, 4, RadType >::handle_non_physical_result<FadType>(const std::string message) const;
325 template FadType PhysicsBase < PHILIP_DIM, 5, RadType >::handle_non_physical_result<FadType>(const std::string message) const;
326 template FadType PhysicsBase < PHILIP_DIM, 6, RadType >::handle_non_physical_result<FadType>(const std::string message) const;
327 template FadType PhysicsBase < PHILIP_DIM, 8, RadType >::handle_non_physical_result<FadType>(const std::string message) const;
328 
329 template FadType PhysicsBase < PHILIP_DIM, 1, FadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
330 template FadType PhysicsBase < PHILIP_DIM, 2, FadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
331 template FadType PhysicsBase < PHILIP_DIM, 3, FadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
332 template FadType PhysicsBase < PHILIP_DIM, 4, FadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
333 template FadType PhysicsBase < PHILIP_DIM, 5, FadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
334 template FadType PhysicsBase < PHILIP_DIM, 6, FadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
335 template FadType PhysicsBase < PHILIP_DIM, 8, FadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
336 
337 template FadType PhysicsBase < PHILIP_DIM, 1, RadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
338 template FadType PhysicsBase < PHILIP_DIM, 2, RadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
339 template FadType PhysicsBase < PHILIP_DIM, 3, RadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
340 template FadType PhysicsBase < PHILIP_DIM, 4, RadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
341 template FadType PhysicsBase < PHILIP_DIM, 5, RadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
342 template FadType PhysicsBase < PHILIP_DIM, 6, RadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
343 template FadType PhysicsBase < PHILIP_DIM, 8, RadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
344 
345 } // Physics namespace
346 } // PHiLiP namespace
347 
PHiLiP::FadFadType
Sacado::Fad::DFad< FadType > FadFadType
Sacado AD type that allows 2nd derivatives.
Definition: ADTypes.hpp:12
PHiLiP::FadType
Sacado::Fad::DFad< double > FadType
Sacado AD type for first derivatives.
Definition: ADTypes.hpp:11
PHiLiP::Physics::PhysicsBase
Base class from which Advection, Diffusion, ConvectionDiffusion, and Euler is derived.
Definition: physics.h:34
PHiLiP::ManufacturedSolutionFunction
Manufactured solution used for grid studies to check convergence orders.
Definition: manufactured_solution.h:22
PHiLiP::RadType
codi_JacobianComputationType RadType
CoDiPaco reverse-AD type for first derivatives.
Definition: ADTypes.hpp:27
PHiLiP::Physics::PhysicsBase::PhysicsBase
PhysicsBase(const Parameters::AllParameters *const parameters_input, const bool has_nonzero_diffusion_input, const bool has_nonzero_physical_source_input, const dealii::Tensor< 2, 3, double > input_diffusion_tensor=Parameters::ManufacturedSolutionParam::get_default_diffusion_tensor(), std::shared_ptr< ManufacturedSolutionFunction< dim, real > > manufactured_solution_function_input=nullptr)
Default constructor that will set the constants.
Definition: physics.cpp:16
PHiLiP
Files for the baseline physics.
Definition: ADTypes.hpp:10
PHiLiP::Parameters::AllParameters
Main parameter class that contains the various other sub-parameter classes.
Definition: all_parameters.h:33
PHiLiP::RadFadType
codi_HessianComputationType RadFadType
Nested reverse-forward mode type for Jacobian and Hessian computation using TapeHelper.
Definition: ADTypes.hpp:28