target_wall_pressure.hpp

#ifndef __PHILIP_TARGET_WALL_PRESSURE_H__
#define __PHILIP_TARGET_WALL_PRESSURE_H__

/* includes */
#include <deal.II/differentiation/ad/sacado_math.h>
#include <deal.II/differentiation/ad/sacado_number_types.h>
#include <deal.II/differentiation/ad/sacado_product_types.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/lac/la_parallel_vector.h>

#include <Sacado.hpp>
#include <iostream>
#include <vector>

#include "dg/dg_base.hpp"
#include "physics/physics.h"
#include "target_functional.h"

namespace PHiLiP {

template <int dim, int nstate, typename real>
class TargetWallPressure : public TargetFunctional<dim, nstate, real>
{
private:
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 


    using TargetFunctional<dim,nstate,real>::evaluate_volume_integrand;
    using TargetFunctional<dim,nstate,real>::evaluate_boundary_integrand;
    using TargetFunctional<dim,nstate,real>::TargetFunctional;

public:

    real evaluate_functional( const bool compute_dIdW = false, const bool compute_dIdX = false, const bool compute_d2I = false) override
    {
        double value = TargetFunctional<dim,nstate,real>::evaluate_functional( compute_dIdW, compute_dIdX, compute_d2I);

        this->pcout << "Target pressure error l2_norm: " << value << "\n";

        return value;
    }



    template<typename real2>
    real2 evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,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<real,nstate> &target_soln_at_q,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &/*soln_grad_at_q*/,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &/*target_soln_grad_at_q*/) const
    {
        if (boundary_id == 1001) {
            assert(soln_at_q.size() == dim+2);
            const Physics::Euler<dim,dim+2,real2> &euler = dynamic_cast< const Physics::Euler<dim,dim+2,real2> &> (physics);

            real2 pressure = euler.compute_pressure (soln_at_q);

            std::array<real2,nstate> target_soln_at_q_real2;
            for (int s=0; s<nstate; ++s) {
                target_soln_at_q_real2[s] = target_soln_at_q[s];
            }
            real2 target_pressure = euler.compute_pressure (target_soln_at_q_real2);
            real2 diff = pressure - target_pressure;
            real2 diff2 = diff*diff;

            return diff2;
        } 
        return (real2) 0.0;
    }


    virtual real evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,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<real,nstate> &target_soln_at_q,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_at_q,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &target_soln_grad_at_q) const override
    {
        return evaluate_boundary_integrand<real>(
            physics,
            boundary_id,
            phys_coord,
            normal,
            target_soln_at_q,
            soln_at_q,
            soln_grad_at_q,
            target_soln_grad_at_q);
    }

    virtual FadFadType evaluate_boundary_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,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<real,nstate> &target_soln_at_q,
        const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &soln_grad_at_q,
        const std::array<dealii::Tensor<1,dim,FadFadType>,nstate> &target_soln_grad_at_q) const override
    {
        return evaluate_boundary_integrand<FadFadType>(
            physics,
            boundary_id,
            phys_coord,
            normal,
            soln_at_q,
            target_soln_at_q,
            soln_grad_at_q,
            target_soln_grad_at_q);
    }


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

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


}; // TargetWallPressure class

} // PHiLiP namespace
#endif