physics.cpp

#include <assert.h>
#include <cmath>
#include <vector>
#include <stdlib.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/mpi.h>
#include <boost/preprocessor/seq/for_each.hpp>

#include "ADTypes.hpp"

#include "physics.h"

namespace PHiLiP {
namespace Physics {

template <int dim, int nspecies, int nstate, typename real>
PhysicsBase<dim,nspecies,nstate,real>::PhysicsBase(
    const Parameters::AllParameters *const                    parameters_input,
    const bool                                                has_nonzero_diffusion_input,
    const bool                                                has_nonzero_physical_source_input,
    const dealii::Tensor<2,3,double>                          input_diffusion_tensor,
    std::shared_ptr< ManufacturedSolutionFunction<dim,nspecies,real> > manufactured_solution_function_input)
    : has_nonzero_diffusion(has_nonzero_diffusion_input)
    , has_nonzero_physical_source(has_nonzero_physical_source_input)
    , all_parameters(parameters_input)
    , non_physical_behavior_type(all_parameters->non_physical_behavior_type)
    , manufactured_solution_function(manufactured_solution_function_input)
    , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0)
{
    // if provided with a null ptr, give it the default manufactured solution
    // currently only necessary for the unit test
    if(!manufactured_solution_function)
            manufactured_solution_function = std::make_shared<ManufacturedSolutionSine<dim,nspecies,real>>(nstate);

    // anisotropic diffusion matrix
    diffusion_tensor[0][0] = input_diffusion_tensor[0][0];
    if constexpr(dim >= 2) {
        diffusion_tensor[0][1] = input_diffusion_tensor[0][1];
        diffusion_tensor[1][0] = input_diffusion_tensor[1][0];
        diffusion_tensor[1][1] = input_diffusion_tensor[1][1];
    }
    if constexpr(dim >= 3) {
        diffusion_tensor[0][2] = input_diffusion_tensor[0][2];
        diffusion_tensor[2][0] = input_diffusion_tensor[2][0];
        diffusion_tensor[1][2] = input_diffusion_tensor[1][2];
        diffusion_tensor[2][1] = input_diffusion_tensor[2][1];
        diffusion_tensor[2][2] = input_diffusion_tensor[2][2];
    }
}

template <int dim, int nspecies, int nstate, typename real>
PhysicsBase<dim,nspecies,nstate,real>::PhysicsBase(
    const Parameters::AllParameters *const                    parameters_input,
    const bool                                                has_nonzero_diffusion_input,
    const bool                                                has_nonzero_physical_source_input,
    std::shared_ptr< ManufacturedSolutionFunction<dim,nspecies,real> > manufactured_solution_function_input)
    : PhysicsBase<dim,nspecies,nstate,real>(
        parameters_input,
        has_nonzero_diffusion_input,
        has_nonzero_physical_source_input,
        Parameters::ManufacturedSolutionParam::get_default_diffusion_tensor(),
        manufactured_solution_function_input)
{ }

template <int dim, int nspecies, int nstate, typename real>
std::array<dealii::Tensor<1,dim,real>,nstate> PhysicsBase<dim,nspecies,nstate,real>::convective_numerical_split_flux (
    const std::array<real,nstate> &/*conservative_soln1*/,
    const std::array<real,nstate> &/*conservative_soln2*/) const
{
    pcout << "ERROR: convective_numerical_split_flux() has not yet been implemented for (overridden by) the selected PDE. Aborting..." <<std::flush;
    std::abort();
    std::array<dealii::Tensor<1,dim,real>,nstate> dummy;
    return dummy;
}

template <int dim, int nspecies, int nstate, typename real>
real PhysicsBase<dim,nspecies,nstate,real>
::max_convective_normal_eigenvalue (
    const std::array<real,nstate> &conservative_soln,
    const dealii::Tensor<1,dim,real> &/*normal*/) const
{
    return max_convective_eigenvalue(conservative_soln);
}

/*
template <int dim, int nspecies, int nstate, typename real>
std::array<dealii::Tensor<1,dim,real>,nstate> PhysicsBase<dim,nspecies,nstate,real>
::artificial_dissipative_flux (
    const real viscosity_coefficient,
    const std::array<real,nstate> &,//solution,
    const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient) const
{
    std::array<dealii::Tensor<1,dim,real>,nstate> diss_flux;
    for (int i=0; i<nstate; i++) {
        for (int d=0; d<dim; d++) {
            diss_flux[i][d] = -viscosity_coefficient*(solution_gradient[i][d]);
        }
    }
    return diss_flux;
}
*/

template <int dim, int nspecies, int nstate, typename real>
std::array<real,nstate> PhysicsBase<dim,nspecies,nstate,real>
::dissipative_flux_dot_normal (
        const std::array<real,nstate> &solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
        const std::array<real,nstate> &filtered_solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &filtered_solution_gradient,
        const bool /*on_boundary*/,
        const dealii::types::global_dof_index cell_index,
        const dealii::Tensor<1,dim,real> &normal,
        const int /*boundary_type*/)
{
    std::array<dealii::Tensor<1,dim,real>,nstate> dissipative_flux = this->dissipative_flux(solution,solution_gradient,filtered_solution,filtered_solution_gradient,cell_index);
    std::array<real,nstate> dissipative_flux_dot_normal;
    dissipative_flux_dot_normal.fill(0.0); // initialize
    for (int s=0; s<nstate; s++) {
        for (int d=0; d<dim; ++d) {
            dissipative_flux_dot_normal[s] += dissipative_flux[s][d] * normal[d];//compute dot product
        }
    }
    return dissipative_flux_dot_normal;
}

template <int dim, int nspecies, int nstate, typename real>
std::array<dealii::Tensor<1,dim,real>,nstate> PhysicsBase<dim,nspecies,nstate,real>
::dissipative_flux (
        const std::array<real,nstate> &solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
        const std::array<real,nstate> &/*filtered_solution*/,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &/*filtered_solution_gradient*/,
        const dealii::types::global_dof_index cell_index)
{
    return this->dissipative_flux(solution,solution_gradient,cell_index);
}

template <int dim, int nspecies, int nstate, typename real>
std::array<real,nstate> PhysicsBase<dim,nspecies,nstate,real>
::artificial_source_term (
    const real viscosity_coefficient,
    const dealii::Point<dim,real> &pos,
    const std::array<real,nstate> &/*solution*/) const
{
    std::array<real,nstate> source;
    
    dealii::Tensor<2,dim,double> artificial_diffusion_tensor;
    for (int i=0;i<dim;i++)
        for (int j=0;j<dim;j++)
            artificial_diffusion_tensor[i][j] = (i==j) ? 1.0 : 0.0;

    for (int istate=0; istate<nstate; istate++) {
        dealii::SymmetricTensor<2,dim,real> manufactured_hessian = this->manufactured_solution_function->hessian (pos, istate);
        //source[istate] = -viscosity_coefficient*scalar_product(artificial_diffusion_tensor,manufactured_hessian);
        source[istate] = 0.0;
        for (int dr=0; dr<dim; ++dr) {
            for (int dc=0; dc<dim; ++dc) {
                source[istate] += artificial_diffusion_tensor[dr][dc] * manufactured_hessian[dr][dc];
            }
        }
        source[istate] *= -viscosity_coefficient;
    }
    return source;
}

template <int dim, int nspecies, int nstate, typename real>
real PhysicsBase<dim,nspecies,nstate,real>
::compute_pressure ( const std::array<real,nstate> &/*conservative_soln*/ ) const
{
    std::cout << "The compute_pressure function has not been implemented for this PDE...Aborting." << std::endl;
    std::abort();
    return 0;
}

template <int dim, int nspecies, int nstate, typename real>
real PhysicsBase<dim,nspecies,nstate,real>
::compute_entropy ( const std::array<real,nstate> &/*conservative_soln*/ ) const
{
    std::cout << "The compute_entropy function has not been implemented for this PDE...Aborting." << std::endl;
    std::abort();
    return 0;
}

template <int dim, int nspecies, int nstate, typename real>
real PhysicsBase<dim,nspecies,nstate,real>
::compute_gamma ( const std::array<real,nstate> &/*conservative_soln*/ ) const
{
    std::cout << "The compute_gamma function has not been implemented for this PDE...Aborting." << std::endl;
    std::abort();
    return 0;
}

template <int dim, int nspecies, int nstate, typename real>
std::array<real,nstate> PhysicsBase<dim,nspecies,nstate,real>
::compute_kinetic_energy_variables ( const std::array<real,nstate> &/*conservative_soln*/ ) const
{
    std::cout << "The compute_kinetic_energy_variables function has not been implemented for this PDE...Aborting." << std::endl;
    std::abort();

    std::array<real,nstate> kinetic_energy_var;
    std::fill(kinetic_energy_var.begin(), kinetic_energy_var.end(), 0.0);

    return kinetic_energy_var;
}

template <int dim, int nspecies, int nstate, typename real>
void PhysicsBase<dim,nspecies,nstate,real>
::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,
        const std::array<real,nstate> &/*filtered_soln_int*/,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &/*filtered_soln_grad_int*/,
        std::array<real,nstate> &soln_bc,
        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
{
    this->boundary_face_values(boundary_type,
                               pos,
                               normal,
                               soln_int,
                               soln_grad_int,
                               soln_bc,
                               soln_grad_bc);
}

template <int dim, int nspecies, int nstate, typename real>
void PhysicsBase<dim,nspecies,nstate,real>
::boundary_face_values_viscous_flux (
        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,
        const std::array<real,nstate> &/*filtered_soln_int*/,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &/*filtered_soln_grad_int*/,
        std::array<real,nstate> &soln_bc,
        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
{
    this->boundary_face_values(boundary_type,
                               pos,
                               normal,
                               soln_int,
                               soln_grad_int,
                               soln_bc,
                               soln_grad_bc);
}

template <int dim, int nspecies, int nstate, typename real>
std::array<real,nstate> PhysicsBase<dim,nspecies,nstate,real>
::physical_source_term (
    const dealii::Point<dim,real> &/*pos*/,
    const std::array<real,nstate> &/*solution*/,
    const std::array<dealii::Tensor<1,dim,real>,nstate> &/*solution_gradient*/,
    const dealii::types::global_dof_index /*cell_index*/) const
{
    std::array<real,nstate> physical_source;
    for (int i=0; i<nstate; i++) {
        physical_source[i] = 0;
    }
    return physical_source;
}

template <int dim, int nspecies, int nstate, typename real>
dealii::Vector<double> PhysicsBase<dim,nspecies,nstate,real>::post_compute_derived_quantities_vector (
    const dealii::Vector<double>              &uh,
    const std::vector<dealii::Tensor<1,dim> > &/*duh*/,
    const std::vector<dealii::Tensor<2,dim> > &/*dduh*/,
    const dealii::Tensor<1,dim>               &/*normals*/,
    const dealii::Point<dim>                  &/*evaluation_points*/) const
{
    dealii::Vector<double> computed_quantities(nstate);
    for (unsigned int s=0; s<nstate; ++s) {
        computed_quantities(s) = uh(s);
    }
    return computed_quantities;
}

template <int dim, int nspecies, int nstate, typename real>
dealii::Vector<double> PhysicsBase<dim,nspecies,nstate,real>::post_compute_derived_quantities_scalar (
    const double              &uh,
    const dealii::Tensor<1,dim> &/*duh*/,
    const dealii::Tensor<2,dim> &/*dduh*/,
    const dealii::Tensor<1,dim> &/*normals*/,
    const dealii::Point<dim>    &/*evaluation_points*/) const
{
    assert(nstate == 1);
    dealii::Vector<double> computed_quantities(nstate);
    for (unsigned int s=0; s<nstate; ++s) {
        computed_quantities(s) = uh;
    }
    return computed_quantities;
}

template <int dim, int nspecies, int nstate, typename real>
std::vector<std::string> PhysicsBase<dim,nspecies,nstate,real> ::post_get_names () const
{
    std::vector<std::string> names;
    for (unsigned int s=0; s<nstate; ++s) {
        std::string varname = "state" + dealii::Utilities::int_to_string(s,1);
        names.push_back(varname);
    }
    return names;
}

template <int dim, int nspecies, int nstate, typename real>
std::vector<dealii::DataComponentInterpretation::DataComponentInterpretation> PhysicsBase<dim,nspecies,nstate,real>
::post_get_data_component_interpretation () const
{
    namespace DCI = dealii::DataComponentInterpretation;
    std::vector<DCI::DataComponentInterpretation> interpretation;
    for (unsigned int s=0; s<nstate; ++s) {
        interpretation.push_back (DCI::component_is_scalar);
    }
    return interpretation;
}

template <int dim, int nspecies, int nstate, typename real>
dealii::UpdateFlags PhysicsBase<dim,nspecies,nstate,real>
::post_get_needed_update_flags () const
{
    return dealii::update_values;
}

template <int dim, int nspecies, int nstate, typename real>
template<typename real2>
real2 PhysicsBase<dim,nspecies,nstate,real>
::handle_non_physical_result(const std::string message) const
{
    if (this->non_physical_behavior_type == NonPhysicalBehaviorEnum::abort_run) {
        std::cout << "ERROR: Non-physical result has been detected. ";
        if (!message.empty()) {
            std::cout << std::endl << "    Message: " << message << std::endl;
        }
        std::cout << " Aborting... " << std::endl << std::flush;
        std::abort();
    } else if (this->non_physical_behavior_type == NonPhysicalBehaviorEnum::print_warning) {
        std::cout << "WARNING: Non-physical result has been detected at a node." << std::endl;
        if (!message.empty()) {
            std::cout << std::endl << "    Message: " << message << std::endl;
        }
    } else if (this->non_physical_behavior_type == NonPhysicalBehaviorEnum::return_big_number) {
        // do nothing -- assume that the test or iterative solver can handle this.
    }
        
    return (real2)BIG_NUMBER;
}
#if PHILIP_SPECIES==1
    // Define a sequence of indices representing the range of nstate
    #define POSSIBLE_NSTATE (1)(2)(3)(4)(5)(6)(8)

    // Define a macro to instantiate functions for a specific nstate
    #define INSTANTIATE_FOR_NSTATE(r, data, nstate) \
        template class PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, double >; \
        template class PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, FadType >; \
        template class PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, RadType >; \
        template class PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType >; \
        template class PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType >; \
        /* -- handle_non_physical_result */ \
        template double PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, double >::handle_non_physical_result<double>(const std::string message) const; \
        template FadType PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, FadType >::handle_non_physical_result<FadType>(const std::string message) const; \
        template RadType PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, RadType >::handle_non_physical_result<RadType>(const std::string message) const; \
        template FadFadType PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType >::handle_non_physical_result<FadFadType>(const std::string message) const; \
        template RadFadType PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType >::handle_non_physical_result<RadFadType>(const std::string message) const; \
        /* instantiate all the real types with real2 = FadType for automatic differentiation in NavierStokes::dissipative_flux_directional_jacobian() */ \
        template FadType PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, double >::handle_non_physical_result<FadType>(const std::string message) const; \
        template FadType PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, RadType >::handle_non_physical_result<FadType>(const std::string message) const; \
        template FadType PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType >::handle_non_physical_result<FadType>(const std::string message) const; \
        template FadType PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
    BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_FOR_NSTATE, _, POSSIBLE_NSTATE)
#else
    #define POSSIBLE_TYPE (double)(FadType)(RadType)(FadFadType)(RadFadType)
    #define INSTANTIATE_TYPES(r, data, type) \
        template class PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type >; \
        template type PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type >::handle_non_physical_result<type>(const std::string message) const;
    BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_TYPES, _, POSSIBLE_TYPE)
    template FadType PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, double >::handle_non_physical_result<FadType>(const std::string message) const;
    template FadType PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, RadType >::handle_non_physical_result<FadType>(const std::string message) const;
    template FadType PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, FadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
    template FadType PhysicsBase < PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, RadFadType >::handle_non_physical_result<FadType>(const std::string message) const;
#endif
} // Physics namespace
} // PHiLiP namespace