RANS fixed bug, renamed file and added boundary conditions

src/dg/strong_dg.cpp

@@ -211,8 +211,8 @@ void DGStrong<dim,nstate,real,MeshType>::assemble_volume_term_derivatives(
 
     std::vector< ADArrayTensor1 > conv_phys_flux_at_q(n_quad_pts);
     std::vector< ADArrayTensor1 > diss_phys_flux_at_q(n_quad_pts);
-    std::vector< ADArray > source_at_q(n_quad_pts);
-    std::vector< ADArray > physical_source_at_q(n_quad_pts);
+    std::vector< ADArray > source_at_q;
+    std::vector< ADArray > physical_source_at_q;
 
     // AD variable
     std::vector< FadType > soln_coeff(n_dofs_cell);
@@ -242,12 +242,14 @@ void DGStrong<dim,nstate,real,MeshType>::assemble_volume_term_derivatives(
         diss_phys_flux_at_q[iquad] = this->pde_physics_fad->dissipative_flux (soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
 
         if(this->pde_physics_fad->has_nonzero_physical_source){
+            physical_source_at_q.resize(n_quad_pts);
             const dealii::Point<dim,real> real_quad_points = fe_values_vol.quadrature_point(iquad);
             dealii::Point<dim,FadType> ad_points;
             for (int d=0;d<dim;++d) { ad_points[d] = real_quad_points[d]; }
             physical_source_at_q[iquad] = this->pde_physics_fad->physical_source_term (ad_points, soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
         }
         if(this->all_parameters->manufactured_convergence_study_param.manufactured_solution_param.use_manufactured_source_term) {
+            source_at_q.resize(n_quad_pts);
             const dealii::Point<dim,real> real_quad_point = fe_values_vol.quadrature_point(iquad);
             dealii::Point<dim,FadType> ad_point;
             for (int d=0;d<dim;++d) { ad_point[d] = real_quad_point[d]; }
@@ -550,8 +552,8 @@ void DGStrong<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
 
     std::vector< realArrayTensor1 > conv_phys_flux_at_q(n_quad_pts);
     std::vector< realArrayTensor1 > diss_phys_flux_at_q(n_quad_pts);
-    std::vector< realArray > source_at_q(n_quad_pts);
-    std::vector< realArray > physical_source_at_q(n_quad_pts);
+    std::vector< realArray > source_at_q;
+    std::vector< realArray > physical_source_at_q;
 
     // AD variable
     std::vector< realtype > soln_coeff(n_dofs_cell);
@@ -578,10 +580,12 @@ void DGStrong<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
         // Evaluate physical convective flux and source term
         conv_phys_flux_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->convective_flux (soln_at_q[iquad]);
         diss_phys_flux_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->dissipative_flux (soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
-        if(DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->has_nonzero_physical_source){
+        if(this->pde_physics_double->has_nonzero_physical_source){
+            physical_source_at_q.resize(n_quad_pts);
             physical_source_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->physical_source_term (fe_values_vol.quadrature_point(iquad), soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
         }
         if(this->all_parameters->manufactured_convergence_study_param.manufactured_solution_param.use_manufactured_source_term) {
+            source_at_q.resize(n_quad_pts);
             source_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->source_term (fe_values_vol.quadrature_point(iquad), soln_at_q[iquad], current_cell_index);
         }
     }
@@ -637,7 +641,7 @@ void DGStrong<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
             rhs = rhs + fe_values_vol.shape_grad_component(itest,iquad,istate) * diss_phys_flux_at_q[iquad][istate] * JxW[iquad];
 
             // Physical source
-            if(DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->has_nonzero_physical_source){
+            if(this->pde_physics_double->has_nonzero_physical_source){
                 rhs = rhs + fe_values_vol.shape_value_component(itest,iquad,istate) * physical_source_at_q[iquad][istate] * JxW[iquad];
             }
             // Source

src/dg/weak_dg.cpp

@@ -577,8 +577,8 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
 
     std::vector< ADArrayTensor1 > conv_phys_flux_at_q(n_quad_pts);
     std::vector< ADArrayTensor1 > diss_phys_flux_at_q(n_quad_pts);
-    std::vector< doubleArray > source_at_q(n_quad_pts);
-    std::vector< doubleArray > physical_source_at_q(n_quad_pts);
+    std::vector< doubleArray > source_at_q;
+    std::vector< doubleArray > physical_source_at_q;
 
     std::vector< real > soln_coeff(n_soln_dofs_int);
     for (unsigned int idof = 0; idof < n_soln_dofs_int; ++idof) {
@@ -628,7 +628,8 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
         conv_phys_flux_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->convective_flux (soln_at_q[iquad]);
         diss_phys_flux_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->dissipative_flux (soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
 
-        if(DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->has_nonzero_physical_source){
+        if(this->pde_physics_double->has_nonzero_physical_source){
+            physical_source_at_q.resize(n_quad_pts);
             const dealii::Point<dim,real> points = fe_values_vol.quadrature_point(iquad);
             physical_source_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->physical_source_term (points,soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
         }
@@ -639,6 +640,7 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
             }
         }
         if(this->all_parameters->manufactured_convergence_study_param.manufactured_solution_param.use_manufactured_source_term) {
+            source_at_q.resize(n_quad_pts);
             const dealii::Point<dim,real> point = fe_values_vol.quadrature_point(iquad);
             source_at_q[iquad] = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->source_term (point, soln_at_q[iquad], current_cell_index);
             //std::array<real,nstate> artificial_source_at_q = DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->artificial_source_term (artificial_diss_coeff, point, soln_at_q[iquad]);
@@ -678,7 +680,7 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term_explicit(
             //// Note that for diffusion, the negative is defined in the physics_double
             rhs = rhs + fe_values_vol.shape_grad_component(itest,iquad,istate) * diss_phys_flux_at_q[iquad][istate] * JxW[iquad];
             // Physical source
-            if(DGBaseState<dim,nstate,real,MeshType>::pde_physics_double->has_nonzero_physical_source){
+            if(this->pde_physics_double->has_nonzero_physical_source){
                 rhs = rhs + fe_values_vol.shape_value_component(itest,iquad,istate) * physical_source_at_q[iquad][istate] * JxW[iquad];
             }
             // Source
@@ -3623,8 +3625,8 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term(
 
     std::vector< ArrayTensor > conv_phys_flux_at_q(n_quad_pts);
     std::vector< ArrayTensor > diss_phys_flux_at_q(n_quad_pts);
-    std::vector< Array > source_at_q(n_quad_pts);
-    std::vector< Array > physical_source_at_q(n_quad_pts);
+    std::vector< Array > source_at_q;
+    std::vector< Array > physical_source_at_q;
 
     for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
         for (int istate=0; istate<nstate; istate++) {
@@ -3642,6 +3644,7 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term(
         diss_phys_flux_at_q[iquad] = physics.dissipative_flux (soln_at_q[iquad], soln_grad_at_q[iquad], current_cell_index);
 
         if(physics.has_nonzero_physical_source){
+            physical_source_at_q.resize(n_quad_pts);
             dealii::Point<dim,real2> ad_points;
             for (int d=0;d<dim;++d) { ad_points[d] = 0.0;}
             for (unsigned int idof = 0; idof < n_metric_dofs; ++idof) {
@@ -3661,6 +3664,7 @@ void DGWeak<dim,nstate,real,MeshType>::assemble_volume_term(
         }
 
         if(this->all_parameters->manufactured_convergence_study_param.manufactured_solution_param.use_manufactured_source_term) {
+            source_at_q.resize(n_quad_pts);
             dealii::Point<dim,real2> ad_point;
             for (int d=0;d<dim;++d) { ad_point[d] = 0.0;}
             for (unsigned int idof = 0; idof < n_metric_dofs; ++idof) {

src/physics/CMakeLists.txt

@@ -15,7 +15,7 @@ SET(SOURCE
     model.cpp
     large_eddy_simulation.cpp
     reynolds_averaged_navier_stokes.cpp
-    negative_spalart_allmaras.cpp
+    negative_spalart_allmaras_rans_model.cpp
     model_factory.cpp)
 
 foreach(dim RANGE 1 3)

src/physics/convection_diffusion.h

@@ -62,9 +62,7 @@ class ConvectionDiffusion : public PhysicsBase <dim, nstate, real>
             diffusion_scaling_coeff(input_diffusion_coefficient),
             hasConvection(convection), 
             hasDiffusion(diffusion)
-    {
-        //static_assert(nstate<=5, "Physics::ConvectionDiffusion() should be created with nstate<=5");
-    };
+    {};
 
     /// Destructor
     ~ConvectionDiffusion () {};

src/physics/model.cpp

@@ -44,37 +44,37 @@ template <int dim, int nstate, typename real>
 void ModelBase<dim, nstate, real>
 ::boundary_face_values (
    const int /*boundary_type*/,
-   const dealii::Point<dim, real> &pos,
-   const dealii::Tensor<1,dim,real> &normal_int,
-   const std::array<real,nstate> &soln_int,
-   const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+   const dealii::Point<dim, real> &/*pos*/,
+   const dealii::Tensor<1,dim,real> &/*normal_int*/,
+   const std::array<real,nstate> &/*soln_int*/,
+   const std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_int*/,
    std::array<real,nstate> &soln_bc,
    std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
 {
-    std::array<real,nstate> boundary_values;
-    std::array<dealii::Tensor<1,dim,real>,nstate> boundary_gradients;
-    for (int s=0; s<nstate; s++) {
-        boundary_values[s] = this->manufactured_solution_function->value (pos, s);
-        boundary_gradients[s] = this->manufactured_solution_function->gradient (pos, s);
-    }
+    //std::array<real,nstate> boundary_values;
+    //std::array<dealii::Tensor<1,dim,real>,nstate> boundary_gradients;
+    //for (int s=0; s<nstate; s++) {
+    //    boundary_values[s] = this->manufactured_solution_function->value (pos, s);
+    //    boundary_gradients[s] = this->manufactured_solution_function->gradient (pos, s);
+    //}
 
     for (int istate=0; istate<dim+2; ++istate) {
         soln_bc[istate] = 0.0;
         soln_grad_bc[istate] = 0.0;
     }
-    for (int istate=dim+2; istate<nstate; ++istate) {
-
-        std::array<real,nstate> characteristic_dot_n = convective_eigenvalues(boundary_values, normal_int);
-        const bool inflow = (characteristic_dot_n[istate] <= 0.);
-
-        if (inflow) { // Dirichlet boundary condition
-            soln_bc[istate] = boundary_values[istate];
-            soln_grad_bc[istate] = soln_grad_int[istate];
-        } else { // Neumann boundary condition
-            soln_bc[istate] = soln_int[istate];
-            soln_grad_bc[istate] = soln_grad_int[istate];
-        }
-    }
+    //for (int istate=dim+2; istate<nstate; ++istate) {
+//
+    //    std::array<real,nstate> characteristic_dot_n = convective_eigenvalues(boundary_values, normal_int);
+    //    const bool inflow = (characteristic_dot_n[istate] <= 0.);
+//
+    //    if (inflow) { // Dirichlet boundary condition
+    //        soln_bc[istate] = boundary_values[istate];
+    //        soln_grad_bc[istate] = soln_grad_int[istate];
+    //    } else { // Neumann boundary condition
+    //        soln_bc[istate] = soln_int[istate];
+    //        soln_grad_bc[istate] = soln_grad_int[istate];
+    //    }
+    //}
 }
 //----------------------------------------------------------------
 template <int dim, int nstate, typename real>

src/physics/model.h

@@ -78,8 +78,7 @@ class ModelBase
         std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_bc*/) const;
 
     /// Returns current vector solution to be used by PhysicsPostprocessor to output current solution.
-    /** The implementation in this Model base class simply returns the stored solution.
-     */
+    /** The implementation in this Model base class simply returns the stored solution. */
     virtual dealii::Vector<double> post_compute_derived_quantities_vector (
         const dealii::Vector<double>              &uh,
         const std::vector<dealii::Tensor<1,dim> > &/*duh*/,
@@ -88,13 +87,11 @@ class ModelBase
         const dealii::Point<dim>                  &/*evaluation_points*/) const;
 
     /// Returns DataComponentInterpretation of the solution to be used by PhysicsPostprocessor to output current solution.
-    /** Treats every solution state as an independent scalar.
-     */
+    /** Treats every solution state as an independent scalar. */
     virtual std::vector<dealii::DataComponentInterpretation::DataComponentInterpretation> post_get_data_component_interpretation () const;
 
     /// Returns names of the solution to be used by PhysicsPostprocessor to output current solution.
-    /** The implementation in this Model base class simply returns "state(dim+2+0), state(dim+2+1), etc.".
-     */
+    /** The implementation in this Model base class simply returns "state(dim+2+0), state(dim+2+1), etc.". */
     virtual std::vector<std::string> post_get_names () const;
 
     // Quantities needed to be updated by DG for the model -- accomplished by DGBase update_model_variables()

src/physics/model_factory.cpp

@@ -9,7 +9,7 @@
 #include "manufactured_solution.h"
 #include "large_eddy_simulation.h"
 #include "reynolds_averaged_navier_stokes.h"
-#include "negative_spalart_allmaras.h"
+#include "negative_spalart_allmaras_rans_model.h"
 
 namespace PHiLiP {
 namespace Physics {

src/physics/negative_spalart_allmaras.cpp → src/physics/negative_spalart_allmaras_rans_model.cpp

@@ -6,7 +6,7 @@
 
 #include "model.h"
 #include "reynolds_averaged_navier_stokes.h"
-#include "negative_spalart_allmaras.h"
+#include "negative_spalart_allmaras_rans_model.h"
 
 namespace PHiLiP {
 namespace Physics {
@@ -586,6 +586,67 @@ ::compute_production_dissipation_cross_term (
 }
 //----------------------------------------------------------------
 template <int dim, int nstate, typename real>
+void ReynoldsAveragedNavierStokes_SAneg<dim,nstate,real>
+::boundary_wall (
+   std::array<real,nstate> &soln_bc,
+   std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
+{
+    // Wall boundary condition for nu_tilde (working variable of negative SA model)
+    // nu_tilde = 0
+    for (int istate=dim+2; istate<nstate; ++istate) {
+        soln_bc[istate] = 0.0;
+        soln_grad_bc[istate] = 0.0;
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ReynoldsAveragedNavierStokes_SAneg<dim,nstate,real>
+::boundary_outflow (
+   const std::array<real,nstate> &soln_int,
+   const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+   std::array<real,nstate> &soln_bc,
+   std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
+{
+    for (int istate=dim+2; istate<nstate; ++istate) {
+        soln_bc[istate] = soln_int[istate];
+        soln_grad_bc[istate] = soln_grad_int[istate];
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ReynoldsAveragedNavierStokes_SAneg<dim,nstate,real>
+::boundary_inflow (
+   const std::array<real,nstate> &/*soln_int*/,
+   const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+   std::array<real,nstate> &soln_bc,
+   std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
+{
+    const real density_bc = this->navier_stokes_physics->density_inf;
+    const real dynamic_viscosity_coefficient_bc = this->navier_stokes_physics->viscosity_coefficient_inf;
+    const real kinematic_viscosity_coefficient_bc = dynamic_viscosity_coefficient_bc/density_bc;
+    for (int istate=dim+2; istate<nstate; ++istate) {
+        soln_bc[istate] = density_bc*3.0*kinematic_viscosity_coefficient_bc;
+        soln_grad_bc[istate] = soln_grad_int[istate];
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ReynoldsAveragedNavierStokes_SAneg<dim,nstate,real>
+::boundary_farfield (
+   std::array<real,nstate> &soln_bc) const
+{
+   const real density_bc = this->navier_stokes_physics->density_inf;
+   const real dynamic_viscosity_coefficient_bc = this->navier_stokes_physics->viscosity_coefficient_inf;
+   const real kinematic_viscosity_coefficient_bc = dynamic_viscosity_coefficient_bc/density_bc;
+
+   // Farfield boundary condition for nu_tilde (working variable of negative SA model)
+   // nu_tilde = 3.0*kinematic_viscosity_inf to 5.0*kinematic_viscosity_inf
+   for (int istate=dim+2; istate<nstate; ++istate) {
+      soln_bc[istate] = density_bc*3.0*kinematic_viscosity_coefficient_bc;
+   }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
 dealii::Vector<double> ReynoldsAveragedNavierStokes_SAneg<dim,nstate,real>
 ::post_compute_derived_quantities_vector (
     const dealii::Vector<double>              &uh,

src/physics/negative_spalart_allmaras.h → src/physics/negative_spalart_allmaras_rans_model.h

@@ -1,5 +1,5 @@
-#ifndef __NEGATIVE_SPALART_ALLMARAS__
-#define __NEGATIVE_SPALART_ALLMARAS__
+#ifndef __NEGATIVE_SPALART_ALLMARAS_RANS_MODEL__
+#define __NEGATIVE_SPALART_ALLMARAS_RANS_MODEL__
 
 #include "euler.h"
 #include "navier_stokes.h"
@@ -115,6 +115,35 @@ class ReynoldsAveragedNavierStokes_SAneg : public ReynoldsAveragedNavierStokesBa
     template<typename real2> real2 scale_coefficient(
         const real2 coefficient) const;
 
+    /// Wall boundary condition
+    /** Reference: Steven R. Allmaras. (2012). "Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model."
+     *  eq.(7)
+     */ 
+    void boundary_wall (
+        std::array<real,nstate> &soln_bc,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;
+
+    /// Inflow boundary condition
+    void boundary_outflow (
+        const std::array<real,nstate> &soln_int,
+        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+        std::array<real,nstate> &soln_bc,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;
+
+    /// Inflow boundary condition
+    void boundary_inflow (
+        const std::array<real,nstate> &soln_int,
+        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+        std::array<real,nstate> &soln_bc,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;
+
+    /// Farfield boundary conditions based on freestream values
+    /** Reference: Steven R. Allmaras. (2012). "Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model."
+     *  eq.(8)
+     */ 
+    void boundary_farfield (
+        std::array<real,nstate> &soln_bc) const;
+
 private:
     /// Templated nondimensionalized eddy viscosity for the negative SA model.
     /** Reference: Steven R. Allmaras. (2012). "Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model."

src/physics/physics_model.cpp

@@ -248,6 +248,11 @@ ::boundary_face_values (
         std::array<real,nstate_baseline_physics> baseline_soln_bc;
         std::array<dealii::Tensor<1,dim,real>,nstate_baseline_physics> baseline_soln_grad_bc;
 
+        for (int istate=0; istate<nstate_baseline_physics; istate++) {
+            baseline_soln_bc[istate]      = 0;
+            baseline_soln_grad_bc[istate] = 0;
+        }
+
         physics_baseline->boundary_face_values(
                 boundary_type, pos, normal_int, baseline_soln_int, baseline_soln_grad_int, 
                 baseline_soln_bc, baseline_soln_grad_bc);