AD Strong DG - Implicit Inviscid - Computes Jacobians

src/ADTypes.hpp

@@ -26,6 +26,15 @@ using codi_HessianComputationType  = codi::RealReversePrimalIndexGen< codi::Real
 //using RadFadType = codi_JacobianComputationType; ///< Reverse only mode that only allows Jacobian computation.
 using RadType = codi_JacobianComputationType; ///< CoDiPaco reverse-AD type for first derivatives.
 using RadFadType = codi_HessianComputationType ; ///< Nested reverse-forward mode type for Jacobian and Hessian computation using TapeHelper.
+
+/// Tape helper class derived from CoDiPack 
+template<typename adtype>
+class codi_TapeHelper : public codi::TapeHelper<adtype>
+{
+public:
+    /// Returns reference to codi::TapeHelper<adtype>::inputValues which can be used in PHiLiP to add inputs that have already been registered in the global tape.
+    std::vector<typename adtype::GradientData>& getinputValues() {return this->inputValues;} 
+};
 } // PHiLiP namespace
 
 #endif

src/dg/dg_base_state.cpp

@@ -180,7 +180,7 @@ void DGBaseState<dim, nstate, real, MeshType>::update_model_variables() {
 
 template <int dim, int nstate, typename real, typename MeshType>
 void DGBaseState<dim, nstate, real, MeshType>::set_use_auxiliary_eq() {
-    this->use_auxiliary_eq = pde_physics_double->has_nonzero_diffusion;
+    this->use_auxiliary_eq = (pde_physics_double->has_nonzero_diffusion && !this->all_parameters->use_weak_form) ? true : false;
 }
 
 template <int dim, int nstate, typename real, typename MeshType>
@@ -236,4 +236,4 @@ BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_SHARED, _, POSSIBLE_NSTATE)
     template class DGBaseState <PHILIP_DIM, index, double, dealii::Triangulation<PHILIP_DIM>>;
 BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_TRIA, _, POSSIBLE_NSTATE)
 
-}  // namespace PHiLiP
+}  // namespace PHiLiP

src/dg/dg_base_state.hpp

@@ -113,6 +113,7 @@ class DGBaseState : public DGBase<dim, real, MeshType>
     /** Usually called after setting physics.
      */
     void reset_numerical_fluxes();
+
 }; // end of DGBaseState class
 
 }  // namespace PHiLiP

src/flow_solver/flow_solver_cases/limiter_convergence_tests.cpp

@@ -109,8 +109,8 @@ void LimiterConvergenceTests<dim, nstate>::check_limiter_principle(DGBase<dim, d
     const unsigned int init_grid_degree = dg.max_grid_degree;
     const unsigned int poly_degree = this->all_param.flow_solver_param.poly_degree;
     //Constructor for the operators
-    OPERATOR::basis_functions<dim, 2 * dim, double> soln_basis(1, poly_degree, init_grid_degree);
-    OPERATOR::vol_projection_operator<dim, 2 * dim, double> soln_basis_projection_oper(1, dg.max_degree, init_grid_degree);
+    OPERATOR::basis_functions<dim, 2 * dim> soln_basis(1, poly_degree, init_grid_degree);
+    OPERATOR::vol_projection_operator<dim, 2 * dim> soln_basis_projection_oper(1, dg.max_degree, init_grid_degree);
 
 
     // Build the oneD operator to perform interpolation/projection
@@ -217,4 +217,4 @@ void LimiterConvergenceTests<dim, nstate>::compute_unsteady_data_and_write_to_ta
 #endif
 
 }
-}
+}

src/flow_solver/flow_solver_cases/non_periodic_cube_flow.cpp

@@ -69,8 +69,8 @@ void NonPeriodicCubeFlow<dim, nstate>::check_positivity_density(DGBase<dim, doub
     const unsigned int init_grid_degree = dg.max_grid_degree;
     const unsigned int poly_degree = this->all_param.flow_solver_param.poly_degree;
     //Constructor for the operators
-    OPERATOR::basis_functions<dim, 2 * dim, double> soln_basis(1, poly_degree, init_grid_degree);
-    OPERATOR::vol_projection_operator<dim, 2 * dim, double> soln_basis_projection_oper(1, dg.max_degree, init_grid_degree);
+    OPERATOR::basis_functions<dim, 2 * dim> soln_basis(1, poly_degree, init_grid_degree);
+    OPERATOR::vol_projection_operator<dim, 2 * dim> soln_basis_projection_oper(1, dg.max_degree, init_grid_degree);
 
 
     // Build the oneD operator to perform interpolation/projection
@@ -164,4 +164,4 @@ void NonPeriodicCubeFlow<dim, nstate>::compute_unsteady_data_and_write_to_table(
     template class NonPeriodicCubeFlow <PHILIP_DIM,PHILIP_DIM+2>;
 #endif
 } // FlowSolver namespace
-} // PHiLiP namespace
+} // PHiLiP namespace

src/flow_solver/flow_solver_cases/periodic_entropy_tests.cpp

@@ -84,8 +84,8 @@ double PeriodicEntropyTests<dim, nstate>::compute_integrated_quantities(DGBase<d
     const unsigned int n_quad_pts = n_quad_pts_; 
 
     // Construct the basis functions and mapping shape functions.
-    OPERATOR::basis_functions<dim,2*dim,double> soln_basis(1, poly_degree, grid_degree); 
-    OPERATOR::mapping_shape_functions<dim,2*dim,double> mapping_basis(1, poly_degree, grid_degree);
+    OPERATOR::basis_functions<dim,2*dim> soln_basis(1, poly_degree, grid_degree); 
+    OPERATOR::mapping_shape_functions<dim,2*dim> mapping_basis(1, poly_degree, grid_degree);
     // Build basis function volume operator and gradient operator from 1D finite element for 1 state.
     soln_basis.build_1D_volume_operator(dg.oneD_fe_collection_1state[poly_degree], quad_1D);
     soln_basis.build_1D_gradient_operator(dg.oneD_fe_collection_1state[poly_degree], quad_1D);

src/flow_solver/flow_solver_cases/periodic_turbulence.cpp

@@ -178,12 +178,12 @@ void PeriodicTurbulence<dim, nstate>::output_velocity_field(
     dealii::FE_DGQArbitraryNodes<1,1> equidistant_finite_element(vol_quad_equidistant_1D);
 
     const unsigned int init_grid_degree = dg->high_order_grid->fe_system.tensor_degree();
-    OPERATOR::basis_functions<dim,2*dim,double> soln_basis(1, dg->max_degree, init_grid_degree); 
+    OPERATOR::basis_functions<dim,2*dim> soln_basis(1, dg->max_degree, init_grid_degree); 
     soln_basis.build_1D_volume_operator(dg->oneD_fe_collection_1state[dg->max_degree], vol_quad_equidistant_1D);
     soln_basis.build_1D_gradient_operator(dg->oneD_fe_collection_1state[dg->max_degree], vol_quad_equidistant_1D);
 
     // mapping basis for the equidistant node set because we output the physical coordinates
-    OPERATOR::mapping_shape_functions<dim,2*dim,double> mapping_basis_at_equidistant(1, dg->max_degree, init_grid_degree);
+    OPERATOR::mapping_shape_functions<dim,2*dim> mapping_basis_at_equidistant(1, dg->max_degree, init_grid_degree);
     mapping_basis_at_equidistant.build_1D_shape_functions_at_grid_nodes(dg->high_order_grid->oneD_fe_system, dg->high_order_grid->oneD_grid_nodes);
     mapping_basis_at_equidistant.build_1D_shape_functions_at_flux_nodes(dg->high_order_grid->oneD_fe_system, vol_quad_equidistant_1D, dg->oneD_face_quadrature);
 
@@ -340,8 +340,8 @@ void PeriodicTurbulence<dim, nstate>::compute_and_update_integrated_quantities(D
     const unsigned int grid_degree = dg.high_order_grid->fe_system.tensor_degree();
     const unsigned int poly_degree = dg.max_degree;
     // Construct the basis functions and mapping shape functions.
-    OPERATOR::basis_functions<dim,2*dim,double> soln_basis(1, poly_degree, grid_degree); 
-    OPERATOR::mapping_shape_functions<dim,2*dim,double> mapping_basis(1, poly_degree, grid_degree);
+    OPERATOR::basis_functions<dim,2*dim> soln_basis(1, poly_degree, grid_degree); 
+    OPERATOR::mapping_shape_functions<dim,2*dim> mapping_basis(1, poly_degree, grid_degree);
     // Build basis function volume operator and gradient operator from 1D finite element for 1 state.
     soln_basis.build_1D_volume_operator(dg.oneD_fe_collection_1state[poly_degree], quad_extra_1D);
     soln_basis.build_1D_gradient_operator(dg.oneD_fe_collection_1state[poly_degree], quad_extra_1D);
@@ -560,14 +560,14 @@ double PeriodicTurbulence<dim, nstate>::compute_current_integrated_numerical_ent
     const unsigned int n_quad_pts = dg->volume_quadrature_collection[poly_degree].size();
     const unsigned int n_shape_fns = n_dofs_cell / nstate;
 
-    OPERATOR::vol_projection_operator<dim,2*dim,double> vol_projection(1, poly_degree, dg->max_grid_degree);
+    OPERATOR::vol_projection_operator<dim,2*dim> vol_projection(1, poly_degree, dg->max_grid_degree);
     vol_projection.build_1D_volume_operator(dg->oneD_fe_collection_1state[poly_degree], dg->oneD_quadrature_collection[poly_degree]);
 
     // Construct the basis functions and mapping shape functions.
-    OPERATOR::basis_functions<dim,2*dim,double> soln_basis(1, poly_degree, dg->max_grid_degree); 
+    OPERATOR::basis_functions<dim,2*dim> soln_basis(1, poly_degree, dg->max_grid_degree); 
     soln_basis.build_1D_volume_operator(dg->oneD_fe_collection_1state[poly_degree], dg->oneD_quadrature_collection[poly_degree]);
 
-    OPERATOR::mapping_shape_functions<dim,2*dim,double> mapping_basis(1, poly_degree, dg->max_grid_degree);
+    OPERATOR::mapping_shape_functions<dim,2*dim> mapping_basis(1, poly_degree, dg->max_grid_degree);
     mapping_basis.build_1D_shape_functions_at_grid_nodes(dg->high_order_grid->oneD_fe_system, dg->high_order_grid->oneD_grid_nodes);
     mapping_basis.build_1D_shape_functions_at_flux_nodes(dg->high_order_grid->oneD_fe_system, dg->oneD_quadrature_collection[poly_degree], dg->oneD_face_quadrature);
 

src/limiter/maximum_principle_limiter.cpp

@@ -123,8 +123,8 @@ void MaximumPrincipleLimiter<dim, nstate, real>::limit(
     //modal coefficients.
     const unsigned int init_grid_degree = grid_degree;
     //Constructor for the operators
-    OPERATOR::basis_functions<dim, 2 * dim, real> soln_basis(1, max_degree, init_grid_degree);
-    OPERATOR::vol_projection_operator<dim, 2 * dim, real> soln_basis_projection_oper(1, max_degree, init_grid_degree);
+    OPERATOR::basis_functions<dim, 2 * dim> soln_basis(1, max_degree, init_grid_degree);
+    OPERATOR::vol_projection_operator<dim, 2 * dim> soln_basis_projection_oper(1, max_degree, init_grid_degree);
 
 
     // Build the oneD operator to perform interpolation/projection
@@ -229,4 +229,4 @@ template class MaximumPrincipleLimiter <PHILIP_DIM, 3, double>;
 template class MaximumPrincipleLimiter <PHILIP_DIM, 4, double>;
 template class MaximumPrincipleLimiter <PHILIP_DIM, 5, double>;
 template class MaximumPrincipleLimiter <PHILIP_DIM, 6, double>;
-} // PHiLiP namespace
+} // PHiLiP namespace

src/limiter/positivity_preserving_limiter.cpp

@@ -197,8 +197,8 @@ void PositivityPreservingLimiter<dim, nstate, real>::limit(
     const unsigned int init_grid_degree = grid_degree;
 
     // Constructor for the operators
-    OPERATOR::basis_functions<dim, 2 * dim, real> soln_basis(1, max_degree, init_grid_degree);
-    OPERATOR::vol_projection_operator<dim, 2 * dim, real> soln_basis_projection_oper(1, max_degree, init_grid_degree);
+    OPERATOR::basis_functions<dim, 2 * dim> soln_basis(1, max_degree, init_grid_degree);
+    OPERATOR::vol_projection_operator<dim, 2 * dim> soln_basis_projection_oper(1, max_degree, init_grid_degree);
 
     // Build the oneD operator to perform interpolation/projection
     soln_basis.build_1D_volume_operator(oneD_fe_collection_1state[max_degree], oneD_quadrature_collection[max_degree]);
@@ -333,4 +333,4 @@ template class PositivityPreservingLimiter <PHILIP_DIM, 3, double>;
 template class PositivityPreservingLimiter <PHILIP_DIM, 4, double>;
 template class PositivityPreservingLimiter <PHILIP_DIM, 5, double>;
 template class PositivityPreservingLimiter <PHILIP_DIM, 6, double>;
-} // PHiLiP namespace
+} // PHiLiP namespace

src/limiter/tvb_limiter.cpp

@@ -123,7 +123,7 @@ std::array<real, nstate> TVBLimiter<dim, nstate, real>::get_neighbour_cell_avg(
     const dealii::LinearAlgebra::distributed::Vector<double>&       solution,
     const dealii::hp::FECollection<dim>&                            fe_collection,
     const dealii::hp::QCollection<dim>&                             volume_quadrature_collection,
-    OPERATOR::basis_functions<dim, 2 * dim, real>                         soln_basis,
+    OPERATOR::basis_functions<dim, 2 * dim>                         soln_basis,
     const int                                                       poly_degree,
     const std::vector<dealii::types::global_dof_index>&             neigh_dofs_indices,
     const unsigned int                                              n_dofs_neigh_cell)
@@ -193,8 +193,8 @@ void TVBLimiter<dim, nstate, real>::limit(
     //modal coefficients.
     const unsigned int init_grid_degree = grid_degree;
     //Constructor for the operators
-    OPERATOR::basis_functions<dim, 2 * dim, real> soln_basis(1, max_degree, init_grid_degree);
-    OPERATOR::vol_projection_operator<dim, 2 * dim, real> soln_basis_projection_oper(1, max_degree, init_grid_degree);
+    OPERATOR::basis_functions<dim, 2 * dim> soln_basis(1, max_degree, init_grid_degree);
+    OPERATOR::vol_projection_operator<dim, 2 * dim> soln_basis_projection_oper(1, max_degree, init_grid_degree);
 
 
     //build the oneD operator to perform interpolation/projection
@@ -306,4 +306,4 @@ template class TVBLimiter <PHILIP_DIM, 4, double>;
 template class TVBLimiter <PHILIP_DIM, 5, double>;
 template class TVBLimiter <PHILIP_DIM, 6, double>;
 #endif
-} // PHiLiP namespace
+} // PHiLiP namespace

src/limiter/tvb_limiter.h

@@ -38,7 +38,7 @@ class TVBLimiter : public BoundPreservingLimiterState <dim, nstate, real>
         const dealii::LinearAlgebra::distributed::Vector<double>&       solution,
         const dealii::hp::FECollection<dim>&                            fe_collection,
         const dealii::hp::QCollection<dim>&                             volume_quadrature_collection,
-        OPERATOR::basis_functions<dim, 2 * dim, real>                         soln_basis,
+        OPERATOR::basis_functions<dim, 2 * dim>                         soln_basis,
         const int                                                       poly_degree,
         const std::vector<dealii::types::global_dof_index>&             neigh_dofs_indices,
         const unsigned int                                              n_dofs_neigh_cell);

src/ode_solver/relaxation_runge_kutta/root_finding_rrk_ode_solver.cpp

@@ -239,16 +239,16 @@ real RootFindingRRKODESolver<dim,real,MeshType>::compute_integrated_numerical_en
     const unsigned int n_quad_pts = dg->volume_quadrature_collection[poly_degree].size();
     const unsigned int n_shape_fns = n_dofs_cell / nstate;
 
-    OPERATOR::vol_projection_operator<dim,2*dim,double> vol_projection(1, poly_degree, dg->max_grid_degree);
+    OPERATOR::vol_projection_operator<dim,2*dim> vol_projection(1, poly_degree, dg->max_grid_degree);
     vol_projection.build_1D_volume_operator(dg->oneD_fe_collection_1state[poly_degree], 
                                             dg->oneD_quadrature_collection[poly_degree]);
 
     // Construct the basis functions and mapping shape functions.
-    OPERATOR::basis_functions<dim,2*dim,double> soln_basis(1, poly_degree, dg->max_grid_degree); 
+    OPERATOR::basis_functions<dim,2*dim> soln_basis(1, poly_degree, dg->max_grid_degree); 
     soln_basis.build_1D_volume_operator(dg->oneD_fe_collection_1state[poly_degree], 
                                         dg->oneD_quadrature_collection[poly_degree]);
 
-    OPERATOR::mapping_shape_functions<dim,2*dim,double> mapping_basis(1, poly_degree, dg->max_grid_degree);
+    OPERATOR::mapping_shape_functions<dim,2*dim> mapping_basis(1, poly_degree, dg->max_grid_degree);
     mapping_basis.build_1D_shape_functions_at_grid_nodes(dg->high_order_grid->oneD_fe_system, 
                                                          dg->high_order_grid->oneD_grid_nodes);
     mapping_basis.build_1D_shape_functions_at_flux_nodes(dg->high_order_grid->oneD_fe_system, 

src/ode_solver/relaxation_runge_kutta/runge_kutta_store_entropy.cpp

@@ -50,10 +50,10 @@ dealii::LinearAlgebra::distributed::Vector<double> RKNumEntropy<dim,real,MeshTyp
     const unsigned int n_quad_pts = dg->volume_quadrature_collection[poly_degree].size();
     const unsigned int n_shape_fns = n_dofs_cell / nstate;
     //We have to project the vector of entropy variables because the mass matrix has an interpolation from solution nodes built into it.
-    OPERATOR::vol_projection_operator<dim,2*dim,double> vol_projection(1, poly_degree, dg->max_grid_degree);
+    OPERATOR::vol_projection_operator<dim,2*dim> vol_projection(1, poly_degree, dg->max_grid_degree);
     vol_projection.build_1D_volume_operator(dg->oneD_fe_collection_1state[poly_degree], dg->oneD_quadrature_collection[poly_degree]);
 
-    OPERATOR::basis_functions<dim,2*dim,double> soln_basis(1, poly_degree, dg->max_grid_degree);
+    OPERATOR::basis_functions<dim,2*dim> soln_basis(1, poly_degree, dg->max_grid_degree);
     soln_basis.build_1D_volume_operator(dg->oneD_fe_collection_1state[poly_degree], dg->oneD_quadrature_collection[poly_degree]);
 
     dealii::LinearAlgebra::distributed::Vector<double> entropy_var_hat_global(dg->right_hand_side);