lift_drag.hpp

#ifndef __PHILIP_LIFT_DRAG_H__
#define __PHILIP_LIFT_DRAG_H__

#include "functional.h"
#include "parameters/all_parameters.h"
#include "physics/physics_factory.h"
#include "physics/navier_stokes.h"

namespace PHiLiP {

#if PHILIP_DIM==1
template <int dim, int nspecies, int nstate, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, int nspecies, int nstate, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class LiftDragFunctional : public Functional<dim, nspecies, nstate, real, MeshType>
{
public:
    enum Functional_types { lift, drag };
private:
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 


    using Functional<dim,nspecies,nstate,real,MeshType>::evaluate_volume_integrand;

    const Functional_types functional_type;

    const Physics::Euler<dim,nspecies,dim+2,FadFadType> &euler_fad_fad;
    const double angle_of_attack;
    const dealii::Tensor<2,dim,double> rotation_matrix;
    const dealii::Tensor<1,dim,double> lift_vector;
    const dealii::Tensor<1,dim,double> drag_vector;

    dealii::Tensor<1,dim,double> force_vector;


    const double force_dimensionalization_factor;

    const Parameters::AllParameters *const all_parameters; 

    double initialize_force_dimensionalization_factor();

    dealii::Tensor<2,dim,double> initialize_rotation_matrix(const double angle_of_attack);


    dealii::Tensor<1,dim,double> initialize_lift_vector(const dealii::Tensor<2,dim,double> &rotation_matrix);


    dealii::Tensor<1,dim,double> initialize_drag_vector(const dealii::Tensor<2,dim,double> &rotation_matrix);

public:
    LiftDragFunctional(
        std::shared_ptr<DGBase<dim,nspecies,real,MeshType>> dg_input,
        const Functional_types functional_type);

    real evaluate_functional( const bool compute_dIdW = false, const bool compute_dIdX = false, const bool compute_d2I = false) override;

public:

    template<typename real2>
    real2 evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nspecies,nstate,real2> &/*physics*/,
        const unsigned int boundary_id,
        const dealii::Point<dim,real2> &/*phys_coord*/,
        const dealii::Tensor<1,dim,real2> &normal,
        const std::array<real2,nstate> &soln_at_q,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &soln_grad_at_q) const
    {
        if (boundary_id == 1001 || boundary_id == 1006) {
            assert(soln_at_q.size() == dim+2);

            using PDE_enum = Parameters::AllParameters::PartialDifferentialEquation;
            std::shared_ptr< Physics::NavierStokes<dim,nspecies,dim+2,real2> > navier_stokes_physics = std::dynamic_pointer_cast<Physics::NavierStokes<dim,nspecies,dim+2,real2>> (Physics::PhysicsFactory<dim,nspecies,dim+2,real2>::create_Physics(this->all_parameters, PDE_enum::navier_stokes, nullptr));

            // Compute pressure (same as Euler physics)
            const real2 pressure = navier_stokes_physics->compute_pressure (soln_at_q);
            
            // Initialize
            dealii::Tensor<1,dim,real2> viscous_tensor_times_normal;
            for(int i=0; i<dim; i++){
                viscous_tensor_times_normal[i] = 0;
            }
            // add viscous stress tensor contribution if viscous (i.e. not Euler)
            if(this->all_parameters->pde_type != PDE_enum::euler) {
                // Compute viscous stress tensor
                const dealii::Tensor<2,dim,real2> viscous_stress_tensor = navier_stokes_physics->compute_viscous_stress_tensor_from_conservative_templated(soln_at_q, soln_grad_at_q);
                // std::cout<<"Norm of viscous stress tensor = "<<  viscous_stress_tensor[0][0]<<std::endl;
                for (int i=0;i<dim;i++){
                    for (int j=0;j<dim;j++){
                        viscous_tensor_times_normal[i]+= viscous_stress_tensor[i][j]*normal[j];
                    }
                }
            }

            return force_dimensionalization_factor * (pressure * (normal * force_vector) -  viscous_tensor_times_normal*force_vector);
        } 
        return (real2) 0.0;
    }


    virtual real evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nspecies,nstate,real> &physics,
        const unsigned int boundary_id,
        const dealii::Point<dim,real> &phys_coord,
        const dealii::Tensor<1,dim,real> &normal,
        const std::array<real,nstate> &soln_at_q,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_boundary_integrand<real>(
            physics,
            boundary_id,
            phys_coord,
            normal,
            soln_at_q,
            soln_grad_at_q);
    }


    virtual FadFadType evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nspecies,nstate,FadFadType> &physics,
        const unsigned int boundary_id,
        const dealii::Point<dim,FadFadType> &phys_coord,
        const dealii::Tensor<1,dim,FadFadType> &normal,
        const std::array<FadFadType,nstate> &soln_at_q,
        const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &soln_grad_at_q) const override
    {
        return evaluate_boundary_integrand<FadFadType>(
            physics,
            boundary_id,
            phys_coord,
            normal,
            soln_at_q,
            soln_grad_at_q);
    }


    virtual real evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nspecies,nstate,real> &/*physics*/,
        const dealii::Point<dim,real> &/*phys_coord*/,
        const std::array<real,nstate> &/*soln_at_q*/,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_at_q*/) const
    { return (real) 0.0; }


    virtual FadFadType evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nspecies,nstate,FadFadType> &/*physics*/,
        const dealii::Point<dim,FadFadType> &/*phys_coord*/, const std::array<FadFadType,nstate> &/*soln_at_q*/,
        const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &/*soln_grad_at_q*/) const
    { return (FadFadType) 0.0; }


};

// template <int dim, int nspecies, int nstate, typename real>
// class TargetLiftDragFunctional : public LiftDragFunctional<dim,nspecies,nstate,real>
// {
// private:
//     /// Constructor
//     TargetLiftDragFunctional(
//         std::shared_ptr<DGBase<dim,nspecies,real>> dg_input,
//         const Functional_types functional_type
//         const double target_value = -1e200
//         : LiftDragFunctional(dg_input, functional_type)
//         , target_value(target_value)
//     { }
// 
// 
//     real evaluate_functional(
//         const bool compute_dIdW,
//         const bool compute_dIdX,
//         const bool compute_d2I)
//     {
//         real value = LiftDragFunctional<dim,nspecies,nstate,real>::evaluate_functional(compute_dIdW, compute_dIdX, compute_d2I);
// 
//         return value - target_value
//     }
// 
// };
// 
// template <int dim, int nspecies, int nstate, typename real>
// class QuadraticPenaltyTargetLiftDragFunctional : public TargetLiftDragFunctional<dim,nspecies,nstate,real>
// {
// public:
// 
//     double penalty;
// 
//     /// Constructor
//     QuadraticPenaltyTargetLiftDragFunctional(
//         std::shared_ptr<DGBase<dim,nspecies,real>> dg_input,
//         const Functional_types functional_type
//         const double target_value = -1e200
//         const double penalty = 0
//         : TargetLiftDragFunctional(dg_input, functional_type, target_value)
//         , penalty(penalty)
//     { }
// 
// 
//     real evaluate_functional(
//         const bool compute_dIdW,
//         const bool compute_dIdX,
//         const bool compute_d2I)
//     {
//         real value = TargetLiftDragFunctional<dim,nspecies,nstate,real>::evaluate_functional((compute_dIdW || compute_d2I), (compute_dIdX || compute_d2I), compute_d2I);
// 
//         if (compute_dIdW) {
//             const real scaling = 2.0*value;
//             this->dIdw *= scaling;
//         }
// 
//         if (compute_dIdX) {
//             const real scaling = 2.0*value;
//             this->dIdX *= scaling;
//         }
// 
//         if (compute_d2I) {
//             const real scaling = 2.0*value;
//             this->dIdX *= scaling;
//         }
// 
//         return penalty * value * value;
//     }
// 
// };

} // PHiLiP namespace

#endif