target_functional.h

#ifndef __TARGET_FUNCTIONAL_H__
#define __TARGET_FUNCTIONAL_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 "functional.h"
#include "physics/physics.h"

namespace PHiLiP {


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

    using Functional<dim,nstate,real>::dIdw;
    using Functional<dim,nstate,real>::dIdX;
    using Functional<dim,nstate,real>::current_functional_value;
    using Functional<dim,nstate,real>::d2IdWdW;
    using Functional<dim,nstate,real>::d2IdWdX;
    using Functional<dim,nstate,real>::d2IdXdX;

protected:
    using Functional<dim,nstate,real>::dg;
    using Functional<dim,nstate,real>::physics_fad_fad;

    using Functional<dim,nstate,real>::volume_update_flags; 
    using Functional<dim,nstate,real>::face_update_flags; 
    using Functional<dim,nstate,real>::uses_solution_values; 
    using Functional<dim,nstate,real>::uses_solution_gradient; 


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

    // *  us, but is a typical bug that other people have. This 'using' imports the base class function
    // *  to our derived class even though we don't need it.
    // */
    //using Functional<dim,nstate,real>::evaluate_cell_boundary;


    using Functional<dim,nstate,real>::evaluate_volume_cell_functional;


    using Functional<dim,nstate,real>::evaluate_boundary_cell_functional;
    using Functional<dim,nstate,real>::evaluate_boundary_integrand;

protected:
    const dealii::LinearAlgebra::distributed::Vector<real> target_solution;

public:

    TargetFunctional(
        std::shared_ptr<DGBase<dim,real>> dg_input,
        const bool uses_solution_values = true,
        const bool uses_solution_gradient = true);

    TargetFunctional(
        std::shared_ptr<DGBase<dim,real>> dg_input,
        const dealii::LinearAlgebra::distributed::Vector<real> &target_solution,
        const bool uses_solution_values = true,
        const bool uses_solution_gradient = true);


    TargetFunctional(
        std::shared_ptr<DGBase<dim,real>> dg_input,
        std::shared_ptr< Physics::PhysicsBase<dim,nstate,Sacado::Fad::DFad<Sacado::Fad::DFad<real>>> > _physics_fad_fad,
        const bool uses_solution_values = true,
        const bool uses_solution_gradient = true);


    TargetFunctional(
        std::shared_ptr<DGBase<dim,real>> dg_input,
  const dealii::LinearAlgebra::distributed::Vector<real> &target_solution,
        std::shared_ptr< Physics::PhysicsBase<dim,nstate,Sacado::Fad::DFad<Sacado::Fad::DFad<real>>> > _physics_fad_fad,
        const bool uses_solution_values = true,
        const bool uses_solution_gradient = true);

public:

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

    dealii::LinearAlgebra::distributed::Vector<real> evaluate_dIdw_finiteDifferences(
        DGBase<dim,real> &dg, 
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &physics,
        const double stepsize);

    dealii::LinearAlgebra::distributed::Vector<real> evaluate_dIdX_finiteDifferences(
        DGBase<dim,real> &dg, 
        const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &physics,
        const double stepsize);

private:
    template <typename real2>
    real2 evaluate_volume_cell_functional(
        const Physics::PhysicsBase<dim,nstate,real2> &physics,
        const std::vector< real2 > &soln_coeff,
        const std::vector< real > &target_soln_coeff,
        const dealii::FESystem<dim> &fe_solution,
        const std::vector< real2 > &coords_coeff,
        const dealii::FESystem<dim> &fe_metric,
        const dealii::Quadrature<dim> &volume_quadrature) const;

protected:
    virtual real evaluate_volume_cell_functional(
        const Physics::PhysicsBase<dim,nstate,real> &physics,
        const std::vector< real > &soln_coeff,
        const std::vector< real > &target_soln_coeff,
        const dealii::FESystem<dim> &fe_solution,
        const std::vector< real > &coords_coeff,
        const dealii::FESystem<dim> &fe_metric,
        const dealii::Quadrature<dim> &volume_quadrature) const;
    virtual Sacado::Fad::DFad<Sacado::Fad::DFad<real>> evaluate_volume_cell_functional(
        const Physics::PhysicsBase<dim,nstate,Sacado::Fad::DFad<Sacado::Fad::DFad<real>>> &physics_fad_fad,
        const std::vector< Sacado::Fad::DFad<Sacado::Fad::DFad<real>> > &soln_coeff,
        const std::vector< real > &target_soln_coeff,
        const dealii::FESystem<dim> &fe_solution,
        const std::vector< Sacado::Fad::DFad<Sacado::Fad::DFad<real>> > &coords_coeff,
        const dealii::FESystem<dim> &fe_metric,
        const dealii::Quadrature<dim> &volume_quadrature) const;

private:
    template <typename real2>
    real2 evaluate_boundary_cell_functional(
        const Physics::PhysicsBase<dim,nstate,real2> &physics,
        const unsigned int boundary_id,
        const std::vector< real2 > &soln_coeff,
        const std::vector< real > &target_soln_coeff,
        const dealii::FESystem<dim> &fe_solution,
        const std::vector< real2 > &coords_coeff,
        const dealii::FESystem<dim> &fe_metric,
        const dealii::Quadrature<dim-1> &face_quadrature,
        const unsigned int face_number) const;
protected:
    virtual real evaluate_boundary_cell_functional(
        const Physics::PhysicsBase<dim,nstate,real> &physics,
        const unsigned int boundary_id,
        const std::vector< real > &soln_coeff,
        const std::vector< real > &target_soln_coeff,
        const dealii::FESystem<dim> &fe_solution,
        const std::vector< real > &coords_coeff,
        const dealii::FESystem<dim> &fe_metric,
        const dealii::Quadrature<dim-1> &face_quadrature,
        const unsigned int face_number) const;

    virtual Sacado::Fad::DFad<Sacado::Fad::DFad<real>> evaluate_boundary_cell_functional(
        const Physics::PhysicsBase<dim,nstate,Sacado::Fad::DFad<Sacado::Fad::DFad<real>>> &physics_fad_fad,
        const unsigned int boundary_id,
        const std::vector< Sacado::Fad::DFad<Sacado::Fad::DFad<real>> > &soln_coeff,
        const std::vector< real > &target_soln_coeff,
        const dealii::FESystem<dim> &fe_solution,
        const std::vector< Sacado::Fad::DFad<Sacado::Fad::DFad<real>> > &coords_coeff,
        const dealii::FESystem<dim> &fe_metric,
        const dealii::Quadrature<dim-1> &face_quadrature,
        const unsigned int face_number) const;


    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
    {
        real l2error = 0;

        for (int istate=0; istate<nstate; ++istate) {
            l2error += std::pow(soln_at_q[istate] - target_soln_at_q[istate], 2);
        }

        return l2error;
    }

    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
    {
        FadFadType l2error = 0;
        for (int istate=0; istate<nstate; ++istate) {
            l2error += std::pow(soln_at_q[istate] - target_soln_at_q[istate], 2);
        }
        return l2error;
    }


    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
    {
        real l2error = 0;
        for (int istate=0; istate<nstate; ++istate) {
            l2error += std::pow(soln_at_q[istate] - target_soln_at_q[istate], 2);
        }
        return l2error;
    }

    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
    {
        FadFadType l2error = 0;
        for (int istate=0; istate<nstate; ++istate) {
            l2error += std::pow(soln_at_q[istate] - target_soln_at_q[istate], 2);
        }
        return l2error;
    }

    // /// Virtual function for computation of cell boundary functional term
    // /** Used only in the computation of evaluate_function(). If not overriden returns 0. */
    // virtual real evaluate_cell_boundary(
    //     const PHiLiP::Physics::PhysicsBase<dim,nstate,real> &/*physics*/,
    //     const unsigned int /*boundary_id*/,
    //     const dealii::FEFaceValues<dim,dim> &/*fe_values_boundary*/,
    //     std::vector<real> /*soln_coeff*/,
    //     std::vector<real> /*target_soln_coeff*/)
    // {
    //     return (real) 0.0;
    // }

    // /// Virtual function for Sacado computation of cell boundary functional term and derivatives
    // /** Used only in the computation of evaluate_dIdw(). If not overriden returns 0. */
    // virtual FadFadType evaluate_cell_boundary(
    //     const PHiLiP::Physics::PhysicsBase<dim,nstate,FadFadType> &/*physics*/,
    //     const unsigned int /*boundary_id*/,
    //     const dealii::FEFaceValues<dim,dim> &/*fe_values_boundary*/,
    //     std::vector<FadFadType> /*soln_coeff*/,
    //     std::vector<real> /*target_soln_coeff*/)
    // {
    //     return (FadFadType) 0.0;
    // }

}; // TargetFunctional class
} // PHiLiP namespace

#endif // __TARGET_FUNCTIONAL_H__