From 3247262f080e351b290508fa040619e58a29063b Mon Sep 17 00:00:00 2001
From: Julien Brillon <43619715+jbrillon@users.noreply.github.com>
Date: Thu, 29 Sep 2022 16:19:19 -0400
Subject: [PATCH] restructuring ModelBase boundary_face_values

---
 src/physics/model.cpp                         | 174 +++++++++++++++---
 src/physics/model.h                           |  69 ++++++-
 .../negative_spalart_allmaras_rans_model.h    |   8 +-
 .../reynolds_averaged_navier_stokes.cpp       |  47 -----
 src/physics/reynolds_averaged_navier_stokes.h |  35 +---
 5 files changed, 212 insertions(+), 121 deletions(-)

diff --git a/src/physics/model.cpp b/src/physics/model.cpp
index 089eb30f7..60179d528 100644
--- a/src/physics/model.cpp
+++ b/src/physics/model.cpp
@@ -15,6 +15,8 @@ template <int dim, int nstate, typename real>
 ModelBase<dim, nstate, real>::ModelBase(
     std::shared_ptr< ManufacturedSolutionFunction<dim,real> > manufactured_solution_function_input):
         manufactured_solution_function(manufactured_solution_function_input)
+        , mpi_communicator(MPI_COMM_WORLD)
+        , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0)
 { 
     // if provided with a null ptr, give it the default manufactured solution
     // currently only necessary for the unit test
@@ -41,40 +43,158 @@ ::physical_source_term (
 }
 //----------------------------------------------------------------
 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*/,
+void ModelBase<dim,nstate,real>
+::boundary_manufactured_solution (
+    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
+{
+    // Do nothing for nstate==(dim+2)
+    if constexpr(nstate>(dim+2)) {
+        pcout << "Error: boundary_manufactured_solution() not implemented in class derived from ModelBase with nstate>(dim+2)." << std::endl;
+        pcout << "Aborting..." << std::endl;
+        std::abort();
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ModelBase<dim,nstate,real>
+::boundary_wall (
+   std::array<real,nstate> &/*soln_bc*/,
+   std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_bc*/) const
+{
+    // Do nothing for nstate==(dim+2)
+    if constexpr(nstate>(dim+2)) {
+        pcout << "Error: boundary_wall() not implemented in class derived from ModelBase with nstate>(dim+2)." << std::endl;
+        pcout << "Aborting..." << std::endl;
+        std::abort();
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ModelBase<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
+{
+    // Do nothing for nstate==(dim+2)
+    if constexpr(nstate>(dim+2)) {
+        pcout << "Error: boundary_outflow() not implemented in class derived from ModelBase with nstate>(dim+2)." << std::endl;
+        pcout << "Aborting..." << std::endl;
+        std::abort();
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ModelBase<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
+{
+    // Do nothing for nstate==(dim+2)
+    if constexpr(nstate>(dim+2)) {
+        pcout << "Error: boundary_inflow() not implemented in class derived from ModelBase with nstate>(dim+2)." << std::endl;
+        pcout << "Aborting..." << std::endl;
+        std::abort();
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ModelBase<dim,nstate,real>
+::boundary_farfield (
+   std::array<real,nstate> &/*soln_bc*/) const
+{
+    // Do nothing for nstate==(dim+2)
+    if constexpr(nstate>(dim+2)) {
+        pcout << "Error: boundary_farfield() not implemented in class derived from ModelBase with nstate>(dim+2)." << std::endl;
+        pcout << "Aborting..." << std::endl;
+        std::abort();
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ModelBase<dim,nstate,real>
+::boundary_slip_wall (
    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
+{
+    // Do nothing for nstate==(dim+2)
+    if constexpr(nstate>(dim+2)) {
+        pcout << "Error: boundary_slip_wall() not implemented in class derived from ModelBase with nstate>(dim+2)." << std::endl;
+        pcout << "Aborting..." << std::endl;
+        std::abort();
+    }
+}
+//----------------------------------------------------------------
+template <int dim, int nstate, typename real>
+void ModelBase<dim,nstate,real>
+::boundary_riemann (
+   const dealii::Tensor<1,dim,real> &/*normal_int*/,
+   const std::array<real,nstate> &/*soln_int*/,
+   std::array<real,nstate> &/*soln_bc*/) const
+{
+    // Do nothing for nstate==(dim+2)
+    if constexpr(nstate>(dim+2)) {
+        pcout << "Error: boundary_riemann() not implemented in class derived from ModelBase with nstate>(dim+2)." << std::endl;
+        pcout << "Aborting..." << std::endl;
+        std::abort();
+    }
+}
+//----------------------------------------------------------------
+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,
    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);
-    //}
-
-    for (int istate=0; istate<dim+2; ++istate) {
-        soln_bc[istate] = 0.0;
-        soln_grad_bc[istate] = 0.0;
+    if (boundary_type == 1000) {
+        // Manufactured solution boundary condition
+        boundary_manufactured_solution (pos, normal_int, soln_int, soln_grad_int, soln_bc, soln_grad_bc);
+    } 
+    else if (boundary_type == 1001) {
+        // Wall boundary condition for working variables of RANS turbulence model
+        boundary_wall (soln_bc, soln_grad_bc);
+    } 
+    else if (boundary_type == 1002) {
+        // Outflow boundary condition 
+        boundary_outflow (soln_int, soln_grad_int, soln_bc, soln_grad_bc);
+    } 
+    else if (boundary_type == 1003) {
+        // Inflow boundary condition
+        boundary_inflow (soln_int, soln_grad_int, soln_bc, soln_grad_bc);
+    } 
+    else if (boundary_type == 1004) {
+        // Riemann-based farfield boundary condition
+        boundary_riemann (normal_int, soln_int, soln_bc);
+    } 
+    else if (boundary_type == 1005) {
+        // Simple farfield boundary condition
+        boundary_farfield(soln_bc);
+    } 
+    else if (boundary_type == 1006) {
+        // Slip wall boundary condition
+        boundary_slip_wall (normal_int, soln_int, soln_grad_int, soln_bc, soln_grad_bc);
+    } 
+    else {
+        pcout << "Invalid boundary_type: " << boundary_type << " in ModelBase.cpp" << std::endl;
+        std::abort();
     }
-    //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];
-    //    }
-    //}
+    // Note: this does not get called when nstate==dim+2 since baseline physics takes care of it
 }
 //----------------------------------------------------------------
 template <int dim, int nstate, typename real>
diff --git a/src/physics/model.h b/src/physics/model.h
index 6744a9c82..45e0662c2 100644
--- a/src/physics/model.h
+++ b/src/physics/model.h
@@ -1,6 +1,7 @@
 #ifndef __MODEL__
 #define __MODEL__
 
+#include <deal.II/base/conditional_ostream.h>
 #include <deal.II/base/types.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/lac/la_parallel_vector.templates.h>
@@ -27,6 +28,11 @@ class ModelBase
     /// Manufactured solution function
     std::shared_ptr< ManufacturedSolutionFunction<dim,real> > manufactured_solution_function;
 
+protected:
+    const MPI_Comm mpi_communicator; ///< MPI communicator.
+    dealii::ConditionalOStream pcout; ///< Parallel std::cout that only outputs on mpi_rank==0
+
+public:
     /// Convective flux terms additional to the baseline physics (including convective flux terms in additional PDEs of model)
     virtual std::array<dealii::Tensor<1,dim,real>,nstate> 
     convective_flux (
@@ -67,15 +73,15 @@ class ModelBase
         const std::array<real,nstate> &solution,
         const dealii::types::global_dof_index cell_index) const = 0;
 
-    /// Evaluates boundary values and gradients on the other side of the face.
-    virtual void boundary_face_values (
-        const int /*boundary_type*/,
-        const dealii::Point<dim, real> &/*pos*/,
-        const dealii::Tensor<1,dim,real> &/*normal*/,
-        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;
+    /// Boundary condition handler
+    void boundary_face_values (
+        const int boundary_type,
+        const dealii::Point<dim, real> &pos,
+        const dealii::Tensor<1,dim,real> &normal,
+        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;
 
     /// 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. */
@@ -97,6 +103,51 @@ class ModelBase
     // Quantities needed to be updated by DG for the model -- accomplished by DGBase update_model_variables()
     dealii::LinearAlgebra::distributed::Vector<int> cellwise_poly_degree; ///< Cellwise polynomial degree
     dealii::LinearAlgebra::distributed::Vector<double> cellwise_volume; ////< Cellwise element volume
+
+protected:
+    /// Evaluate the manufactured solution boundary conditions.
+    virtual void boundary_manufactured_solution (
+        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;
+
+    /// Wall boundary condition
+    virtual void boundary_wall (
+        std::array<real,nstate> &soln_bc,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;
+
+    /// Outflow Boundary Condition 
+    virtual 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 conditions 
+    virtual 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
+    virtual void boundary_farfield (
+        std::array<real,nstate> &soln_bc) const;
+
+    virtual void boundary_riemann (
+        const dealii::Tensor<1,dim,real> &normal_int,
+        const std::array<real,nstate> &soln_int,
+        std::array<real,nstate> &soln_bc) const;
+
+    virtual void boundary_slip_wall (
+        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;
 };
 
 } // Physics namespace
diff --git a/src/physics/negative_spalart_allmaras_rans_model.h b/src/physics/negative_spalart_allmaras_rans_model.h
index ee66fe30f..48689b173 100644
--- a/src/physics/negative_spalart_allmaras_rans_model.h
+++ b/src/physics/negative_spalart_allmaras_rans_model.h
@@ -121,28 +121,28 @@ class ReynoldsAveragedNavierStokes_SAneg : public ReynoldsAveragedNavierStokesBa
      */ 
     void boundary_wall (
         std::array<real,nstate> &soln_bc,
-        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const override;
 
     /// 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;
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const override;
 
     /// 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;
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const override;
 
     /// 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;
+        std::array<real,nstate> &soln_bc) const override;
 
 private:
     /// Templated nondimensionalized eddy viscosity for the negative SA model.
diff --git a/src/physics/reynolds_averaged_navier_stokes.cpp b/src/physics/reynolds_averaged_navier_stokes.cpp
index 1c41f4daf..249c5925c 100644
--- a/src/physics/reynolds_averaged_navier_stokes.cpp
+++ b/src/physics/reynolds_averaged_navier_stokes.cpp
@@ -734,53 +734,6 @@ ::boundary_manufactured_solution (
     }
 }
 //----------------------------------------------------------------
-template <int dim, int nstate, typename real>
-void ReynoldsAveragedNavierStokesBase<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,
-   std::array<real,nstate> &soln_bc,
-   std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
-{
-    if (boundary_type == 1000) {
-        // Manufactured solution boundary condition
-        boundary_manufactured_solution (pos, normal_int, soln_int, soln_grad_int, soln_bc, soln_grad_bc);
-    } 
-    else if (boundary_type == 1001) {
-        // Wall boundary condition for working variables of RANS turbulence model
-        boundary_wall (soln_bc, soln_grad_bc);
-    } 
-    else if (boundary_type == 1002) {
-        // Outflow boundary condition 
-        boundary_outflow (soln_int, soln_grad_int, soln_bc, soln_grad_bc);
-    } 
-    else if (boundary_type == 1003) {
-        // Inflow boundary condition
-        boundary_inflow (soln_int, soln_grad_int, soln_bc, soln_grad_bc);
-    } 
-    else if (boundary_type == 1004) {
-        // Riemann-based farfield boundary condition
-        std::cout << "Riemann boundary condition is not implemented!" << std::endl;
-        std::abort();
-    } 
-    else if (boundary_type == 1005) {
-        // Simple farfield boundary condition
-        boundary_farfield(soln_bc);
-    } 
-    else if (boundary_type == 1006) {
-        // Slip wall boundary condition
-        std::cout << "Slip wall boundary condition is not implemented!" << std::endl;
-        std::abort();
-    } 
-    else {
-        std::cout << "Invalid boundary_type: " << boundary_type << std::endl;
-        std::abort();
-    }
-}
-//----------------------------------------------------------------
 //----------------------------------------------------------------
 //----------------------------------------------------------------
 // Instantiate explicitly
diff --git a/src/physics/reynolds_averaged_navier_stokes.h b/src/physics/reynolds_averaged_navier_stokes.h
index 3d82a212b..4a46b1b1a 100644
--- a/src/physics/reynolds_averaged_navier_stokes.h
+++ b/src/physics/reynolds_averaged_navier_stokes.h
@@ -120,16 +120,6 @@ class ReynoldsAveragedNavierStokesBase : public ModelBase <dim, nstate, real>
         const std::array<real,nstate> &conservative_solution,
         const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient) const = 0;
 
-    /// Boundary condition handler
-    void boundary_face_values (
-        const int /*boundary_type*/,
-        const dealii::Point<dim, real> &/*pos*/,
-        const dealii::Tensor<1,dim,real> &/*normal*/,
-        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 override;
-
 protected:
     /// Returns the square of the magnitude of the vector 
     template<typename real2> 
@@ -254,30 +244,7 @@ class ReynoldsAveragedNavierStokesBase : public ModelBase <dim, nstate, real>
         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;
-
-    /// Wall boundary condition
-    virtual void boundary_wall (
-        std::array<real,nstate> &soln_bc,
-        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const = 0;
-
-    /// Outflow Boundary Condition 
-    virtual 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 = 0;
-
-    /// Inflow boundary conditions 
-    virtual 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 = 0;
-
-    /// Farfield boundary conditions based on freestream values
-    virtual void boundary_farfield (
-        std::array<real,nstate> &soln_bc) const = 0;
+        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const override;
 };
 
 } // Physics namespace