Skip to content

Commit

Permalink
restructuring ModelBase boundary_face_values
Browse files Browse the repository at this point in the history
  • Loading branch information
@jbrillon
jbrillon committed Sep 29, 2022
1 parent eca3b0f commit 3247262
Show file tree
Hide file tree
Showing 5 changed files with 212 additions and 121 deletions.
174 changes: 147 additions & 27 deletions src/physics/model.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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>
Expand Down
69 changes: 60 additions & 9 deletions src/physics/model.h
View file
Original file line number Diff line number Diff line change
@@ -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>
Expand All @@ -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 (
Expand Down Expand Up @@ -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. */
Expand All @@ -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
Expand Down
8 changes: 4 additions & 4 deletions src/physics/negative_spalart_allmaras_rans_model.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
47 changes: 0 additions & 47 deletions src/physics/reynolds_averaged_navier_stokes.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit 3247262

Please sign in to comment.