Merge after MPI & hp addition

.gitignore

@@ -2,10 +2,8 @@ CMakeCache.txt
 CMakeFiles
 CMakeScripts
 *.cmake
-
-build*
+build/*
 TAGS
-
 Testing
 Makefile
 cmake_install.cmake

src/dg/dg.h

@@ -61,6 +61,8 @@ class DGBase
      *  DGBase cannot use nstate as a compile-time known.  */
     const unsigned int max_degree;
 
+
+
     /// Pointer to all parameters
     const Parameters::AllParameters *const all_parameters;
 
@@ -249,7 +251,7 @@ class DGBase
     std::tuple< dealii::hp::MappingCollection<dim>, dealii::hp::FECollection<dim>,
                 dealii::hp::QCollection<dim>, dealii::hp::QCollection<dim-1>, dealii::hp::QCollection<1>,
                 dealii::hp::FECollection<dim> >
-                create_collection_tuple(const unsigned int max_degree, const int nstate) const;
+                create_collection_tuple(const unsigned int max_degree, const int nstate, const Parameters::AllParameters *const parameters_input) const;
 
     /// Evaluate the integral over the cell volume
     virtual void assemble_volume_terms_implicit(
@@ -368,6 +370,8 @@ class DGWeak : public DGBase<dim, real>
     ~DGWeak();
 
 private:
+
+
     /// Contains the physics of the PDE
     std::shared_ptr < Physics::PhysicsBase<dim, nstate, Sacado::Fad::DFad<real> > > pde_physics;
     /// Convective numerical flux
@@ -454,14 +458,14 @@ class DGStrong : public DGBase<dim, real>
     /// Dissipative numerical flux
     NumericalFlux::NumericalFluxDissipative<dim, nstate, Sacado::Fad::DFad<real> > *diss_num_flux;
 
+
     /// Contains the physics of the PDE
     std::shared_ptr < Physics::PhysicsBase<dim, nstate, real > > pde_physics_double;
     /// Convective numerical flux
     NumericalFlux::NumericalFluxConvective<dim, nstate, real > *conv_num_flux_double;
     /// Dissipative numerical flux
     NumericalFlux::NumericalFluxDissipative<dim, nstate, real > *diss_num_flux_double;
 
-
     /// Evaluate the integral over the cell volume
     void assemble_volume_terms_implicit(
         const dealii::FEValues<dim,dim> &fe_values_volume,

src/dg/strong_dg.cpp

@@ -50,6 +50,7 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_implicit(
     const std::vector<dealii::types::global_dof_index> &cell_dofs_indices,
     dealii::Vector<real> &local_rhs_int_cell)
 {
+
     using ADtype = Sacado::Fad::DFad<real>;
     using ADArray = std::array<ADtype,nstate>;
     using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,ADtype>, nstate >;
@@ -98,6 +99,7 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_implicit(
         // Evaluate physical convective flux and source term
         conv_phys_flux_at_q[iquad] = pde_physics->convective_flux (soln_at_q[iquad]);
         diss_phys_flux_at_q[iquad] = pde_physics->dissipative_flux (soln_at_q[iquad], soln_grad_at_q[iquad]);
+
         if(this->all_parameters->manufactured_convergence_study_param.use_manufactured_source_term) {
             source_at_q[iquad] = pde_physics->source_term (fe_values_vol.quadrature_point(iquad), soln_at_q[iquad]);
         }
@@ -108,12 +110,17 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_implicit(
     // Since we have nodal values of the flux, we use the Lagrange polynomials to obtain the gradients at the quadrature points.
     const dealii::FEValues<dim,dim> &fe_values_lagrange = this->fe_values_collection_volume_lagrange.get_present_fe_values();
     std::vector<ADArray> flux_divergence(n_quad_pts);
+
+    std::array<std::array<std::vector<ADtype>,nstate>,dim> f;
+    std::array<std::array<std::vector<ADtype>,nstate>,dim> g;
+
     for (int istate = 0; istate<nstate; ++istate) {
         for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
             flux_divergence[iquad][istate] = 0.0;
             for ( unsigned int flux_basis = 0; flux_basis < n_quad_pts; ++flux_basis ) {
                 flux_divergence[iquad][istate] += conv_phys_flux_at_q[flux_basis][istate] * fe_values_lagrange.shape_grad(flux_basis,iquad);
             }
+
         }
     }
 
@@ -129,13 +136,15 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_implicit(
 
         ADtype rhs = 0;
 
+
         const unsigned int istate = fe_values_vol.get_fe().system_to_component_index(itest).first;
 
         for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
 
             // Convective
             // Now minus such 2 integrations by parts
             assert(JxW[iquad] - fe_values_lagrange.JxW(iquad) < 1e-14);
+
             rhs = rhs - fe_values_vol.shape_value_component(itest,iquad,istate) * flux_divergence[iquad][istate] * JxW[iquad];
 
             //// Diffusive
@@ -470,6 +479,7 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_explicit(
     const std::vector<dealii::types::global_dof_index> &cell_dofs_indices,
     dealii::Vector<real> &local_rhs_int_cell)
 {
+	//std::cout << "assembling cell terms" << std::endl;
     using ADtype = Sacado::Fad::DFad<real>;
     using ADArray = std::array<ADtype,nstate>;
     using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,ADtype>, nstate >;
@@ -532,7 +542,14 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_explicit(
         for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
             flux_divergence[iquad][istate] = 0.0;
             for ( unsigned int flux_basis = 0; flux_basis < n_quad_pts; ++flux_basis ) {
-                flux_divergence[iquad][istate] += conv_phys_flux_at_q[flux_basis][istate] * fe_values_lagrange.shape_grad(flux_basis,iquad);
+            	if (this->all_parameters->use_split_form == true)
+            	{
+            		flux_divergence[iquad][istate] += 2* pde_physics->convective_numerical_split_flux(soln_at_q[iquad],soln_at_q[flux_basis])[istate] *  fe_values_lagrange.shape_grad(flux_basis,iquad);
+            	}
+            	else
+            	{
+                    flux_divergence[iquad][istate] += conv_phys_flux_at_q[flux_basis][istate] * fe_values_lagrange.shape_grad(flux_basis,iquad);
+            	}
             }
         }
     }
@@ -556,6 +573,7 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_explicit(
             // Convective
             // Now minus such 2 integrations by parts
             assert(JxW[iquad] - fe_values_lagrange.JxW(iquad) < 1e-14);
+
             rhs = rhs - fe_values_vol.shape_value_component(itest,iquad,istate) * flux_divergence[iquad][istate] * JxW[iquad];
 
             //// Diffusive
@@ -720,6 +738,7 @@ void DGStrong<dim,nstate,real>::assemble_face_term_explicit(
     dealii::Vector<real>          &local_rhs_int_cell,
     dealii::Vector<real>          &local_rhs_ext_cell)
 {
+	//std::cout << "assembling face terms" << std::endl;
     using ADtype = Sacado::Fad::DFad<real>;
     using ADArray = std::array<ADtype,nstate>;
     using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,ADtype>, nstate >;
@@ -808,11 +827,16 @@ void DGStrong<dim,nstate,real>::assemble_face_term_explicit(
         //std::cout << "Energy ext" << soln_ext[iquad][nstate-1] << std::endl;
 
         // Evaluate physical convective flux, physical dissipative flux, and source term
+
+        //std::cout <<"evaluating numerical fluxes" <<std::endl;
         conv_num_flux_dot_n[iquad] = conv_num_flux->evaluate_flux(soln_int[iquad], soln_ext[iquad], normal_int);
 
         conv_phys_flux_int[iquad] = pde_physics->convective_flux (soln_int[iquad]);
         conv_phys_flux_ext[iquad] = pde_physics->convective_flux (soln_ext[iquad]);
 
+       // std::cout <<"done evaluating numerical fluxes" <<std::endl;
+
+
         diss_soln_num_flux[iquad] = diss_num_flux->evaluate_solution_flux(soln_int[iquad], soln_ext[iquad], normal_int);
 
         ADArrayTensor1 diss_soln_jump_int, diss_soln_jump_ext;

src/numerical_flux/CMakeLists.txt

@@ -1,6 +1,7 @@
 set(FLUX_SOURCE
     viscous_numerical_flux.cpp
     numerical_flux.cpp
+    split_form_numerical_flux.cpp
     )
 
 foreach(dim RANGE 1 3)

src/numerical_flux/numerical_flux.cpp

@@ -1,7 +1,8 @@
 #include <Sacado.hpp>
 #include <deal.II/differentiation/ad/sacado_product_types.h>
 #include "numerical_flux.h"
-
+#include "viscous_numerical_flux.h"
+#include "split_form_numerical_flux.h"
 
 namespace PHiLiP {
 namespace NumericalFlux {
@@ -34,6 +35,9 @@ ::create_convective_numerical_flux(
     } else if(conv_num_flux_type == AllParam::roe) {
         if constexpr (dim+2==nstate) return new Roe<dim, nstate, real>(physics_input);
     }
+    else if (conv_num_flux_type == AllParam::split_form) {
+    	return new SplitFormNumFlux<dim, nstate, real>(physics_input);
+    }
 
     std::cout << "Invalid numerical flux" << std::endl;
     return nullptr;

src/numerical_flux/split_form_numerical_flux.cpp

@@ -0,0 +1,49 @@
+#include "split_form_numerical_flux.h"
+#include <Sacado.hpp>
+#include <deal.II/differentiation/ad/sacado_product_types.h>
+
+namespace PHiLiP {
+namespace NumericalFlux {
+
+template <int dim, int nstate, typename real>
+std::array<real, nstate> SplitFormNumFlux<dim,nstate,real>::evaluate_flux(
+	const std::array<real, nstate> &soln_int,
+    const std::array<real, nstate> &soln_ext,
+    const dealii::Tensor<1,dim,real> &normal_int) const
+	{
+		//std::cout << "evaluating the split form flux" <<std::endl;
+		using RealArrayVector = std::array<dealii::Tensor<1,dim,real>,nstate>;
+	    RealArrayVector conv_phys_split_flux;
+
+	    conv_phys_split_flux = pde_physics->convective_numerical_split_flux (soln_int,soln_ext);
+	    //std::cout << "done evaluating the conv num split flux" <<std::endl;
+
+	    const real conv_max_eig_int = pde_physics->max_convective_eigenvalue(soln_int);
+	   // std::cout << "1st eig" << std::endl;
+	    const real conv_max_eig_ext = pde_physics->max_convective_eigenvalue(soln_ext);
+	    //std::cout << "2nd eig" << std::endl;
+	    const real conv_max_eig = std::max(conv_max_eig_int, conv_max_eig_ext);
+
+	   // std::cout << "obtained the max eig" <<std::endl;
+	    // Scalar dissipation
+	    std::array<real, nstate> numerical_flux_dot_n;
+	    for (int s=0; s<nstate; s++) {
+	        numerical_flux_dot_n[s] = conv_phys_split_flux[s]*normal_int - 0.5 * conv_max_eig * (soln_ext[s]-soln_int[s]);
+	    }
+	   // std::cout << "about to return split num flux" <<std::endl;
+	    return numerical_flux_dot_n;
+	}
+
+template class SplitFormNumFlux<PHILIP_DIM, 1, double>;
+template class SplitFormNumFlux<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
+template class SplitFormNumFlux<PHILIP_DIM, 2, double>;
+template class SplitFormNumFlux<PHILIP_DIM, 2, Sacado::Fad::DFad<double> >;
+template class SplitFormNumFlux<PHILIP_DIM, 3, double>;
+template class SplitFormNumFlux<PHILIP_DIM, 3, Sacado::Fad::DFad<double> >;
+template class SplitFormNumFlux<PHILIP_DIM, 4, double>;
+template class SplitFormNumFlux<PHILIP_DIM, 4, Sacado::Fad::DFad<double> >;
+template class SplitFormNumFlux<PHILIP_DIM, 5, double>;
+template class SplitFormNumFlux<PHILIP_DIM, 5, Sacado::Fad::DFad<double> >;
+
+}
+}

src/numerical_flux/split_form_numerical_flux.h

@@ -0,0 +1,40 @@
+#ifndef __SPLIT_NUM_FLUX__
+#define __SPLIT_NUM_FLUX__
+
+#include <deal.II/base/tensor.h>
+#include "physics/physics.h"
+#include "numerical_flux/numerical_flux.h"
+
+namespace PHiLiP {
+namespace NumericalFlux {
+
+/// Lax-Friedrichs numerical flux. Derived from NumericalFluxConvective.
+template<int dim, int nstate, typename real>
+class SplitFormNumFlux: public NumericalFluxConvective<dim, nstate, real>
+{
+public:
+
+/// Constructor
+SplitFormNumFlux(std::shared_ptr <Physics::PhysicsBase<dim, nstate, real>> physics_input)
+:
+pde_physics(physics_input)
+{};
+/// Destructor
+~SplitFormNumFlux() {};
+
+/// Returns the convective numerical flux at an interface.
+std::array<real, nstate> evaluate_flux (
+    const std::array<real, nstate> &soln_int,
+    const std::array<real, nstate> &soln_ext,
+    const dealii::Tensor<1,dim,real> &normal1) const;
+
+protected:
+/// Numerical flux requires physics to evaluate convective eigenvalues.
+const std::shared_ptr < Physics::PhysicsBase<dim, nstate, real> > pde_physics;
+
+};
+
+}
+}
+
+#endif

src/ode_solver/ode_solver.cpp

@@ -96,8 +96,9 @@ int ODESolver<dim,real>::steady_state ()
                   << std::endl;
 
         if ((ode_param.ode_output) == Parameters::OutputEnum::verbose &&
-            (this->current_iteration%ode_param.print_iteration_modulo) == 0 )
-        pcout << " Evaluating right-hand side and setting system_matrix to Jacobian... " << std::endl;
+            (this->current_iteration%ode_param.print_iteration_modulo) == 0 ) {
+		pcout << " Evaluating right-hand side and setting system_matrix to Jacobian... " << std::endl;
+	}
 
         this->dg->assemble_residual ();
         this->residual_norm = this->dg->get_residual_l2norm();
@@ -120,7 +121,6 @@ int ODESolver<dim,real>::steady_state ()
                 this->dg->output_results_vtk(file_number);
             }
         }
-
     }
     return 1;
 }
@@ -146,21 +146,29 @@ int ODESolver<dim,real>::advance_solution_time (double time_advance)
     while (this->current_iteration < number_of_time_steps)
     {
         if ((ode_param.ode_output) == Parameters::OutputEnum::verbose &&
-            (this->current_iteration%ode_param.print_iteration_modulo) == 0 )
-        pcout << " ********************************************************** "
-                  << std::endl
-                  << " Iteration: " << this->current_iteration + 1
-                  << " out of: " << number_of_time_steps
-                  << std::endl;
+            (this->current_iteration%ode_param.print_iteration_modulo) == 0 ) {
+		pcout << " ********************************************************** "
+			  << std::endl
+			  << " Iteration: " << this->current_iteration + 1
+			  << " out of: " << number_of_time_steps
+			  << std::endl;
+	}
+        dg->assemble_residual(false);
 
         if ((ode_param.ode_output) == Parameters::OutputEnum::verbose &&
-            (this->current_iteration%ode_param.print_iteration_modulo) == 0 )
-        pcout << " Evaluating right-hand side and setting system_matrix to Jacobian... " << std::endl;
+            (this->current_iteration%ode_param.print_iteration_modulo) == 0 ) {
+		pcout << " Evaluating right-hand side and setting system_matrix to Jacobian... " << std::endl;
+	}
 
         step_in_time(constant_time_step);
 
+
+   if (this->current_iteration%ode_param.print_iteration_modulo == 0) {
+		this->dg->output_results_vtk(this->current_iteration);
+	}
         ++(this->current_iteration);
 
+
         //this->dg->output_results_vtk(this->current_iteration);
     }
     return 1;
@@ -181,8 +189,9 @@ void Implicit_ODESolver<dim,real>::step_in_time (real dt)
     this->dg->add_mass_matrices(1.0/dt);
 
     if ((ode_param.ode_output) == Parameters::OutputEnum::verbose &&
-        (this->current_iteration%ode_param.print_iteration_modulo) == 0 )
-    pcout << " Evaluating system update... " << std::endl;
+        (this->current_iteration%ode_param.print_iteration_modulo) == 0 ) {
+	    pcout << " Evaluating system update... " << std::endl;
+    }
 
     solve_linear (
         this->dg->system_matrix,
@@ -200,7 +209,7 @@ void Explicit_ODESolver<dim,real>::step_in_time (real dt)
 {
     // this->dg->assemble_residual (); // Not needed since it is called in the base class for time step
     this->current_time += dt;
-    const int rk_order = 3;
+    const int rk_order = 1;
     if (rk_order == 1) {
         this->dg->global_inverse_mass_matrix.vmult(this->solution_update, this->dg->right_hand_side);
         this->update_norm = this->solution_update.l2_norm();