diffusion_exact_adjoint.h

#ifndef __DIFFUSION_EXACT_ADJOINT_H__
#define __DIFFUSION_EXACT_ADJOINT_H__

#include <deal.II/base/tensor.h>
#include <deal.II/base/types.h>

#include <memory>

#include "dg/dg_base.hpp"
#include "functional/functional.h"
#include "parameters/all_parameters.h"
#include "physics/convection_diffusion.h"
#include "physics/manufactured_solution.h"
#include "physics/physics.h"
#include "tests.h"

namespace PHiLiP {
namespace Tests {

template <int dim, int nspecies, typename real>
class ManufacturedSolutionU : public ManufacturedSolutionFunction <dim, nspecies, real>
{
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    real value (const dealii::Point<dim,real> &pos, const unsigned int istate = 0) const override;

    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &pos, const unsigned int istate = 0) const override;

    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &/* point */, const unsigned int /* istate = 0 */) const 
    {
        return dealii::SymmetricTensor<2,dim,real>();
    }
};

template <int dim, int nspecies, typename real>
class ManufacturedSolutionV : public ManufacturedSolutionFunction <dim, nspecies, real>
{
protected:
    using dealii::Function<dim,real>::value;
    using dealii::Function<dim,real>::gradient;
    using dealii::Function<dim,real>::hessian;
public:
    real value (const dealii::Point<dim,real> &pos, const unsigned int istate = 0) const override;

    dealii::Tensor<1,dim,real> gradient (const dealii::Point<dim,real> &pos, const unsigned int istate = 0) const override;

    dealii::SymmetricTensor<2,dim,real> hessian (const dealii::Point<dim,real> &/* point */, const unsigned int /* istate = 0 */) const override
    {
        return dealii::SymmetricTensor<2,dim,real>();
    }
};

template <int dim, int nspecies, int nstate, typename real>
class diffusion_objective : public Physics::ConvectionDiffusion <dim, nspecies, nstate, real>
{
protected:
    // For overloading the virtual functions defined in PhysicsBase
    using Physics::PhysicsBase<dim,nspecies,nstate,real>::dissipative_flux;
    using Physics::PhysicsBase<dim,nspecies,nstate,real>::source_term;
public:
    diffusion_objective(
        const Parameters::AllParameters *const                  parameters_input,
        const bool                                              convection, 
        const bool                                              diffusion,
        std::shared_ptr<ManufacturedSolutionFunction<dim,nspecies,real>> manufactured_solution_function): 
            Physics::ConvectionDiffusion<dim,nspecies,nstate,real>::ConvectionDiffusion(
                parameters_input,
                convection, 
                diffusion,
                default_diffusion_tensor(),
                Parameters::ManufacturedSolutionParam::get_default_advection_vector(),
                default_diffusion_coefficient(),
                manufactured_solution_function)
    {}

    static double default_diffusion_coefficient()
    {
        return -1.0;
    }

    static dealii::Tensor<2,3,double> default_diffusion_tensor()
    {
        dealii::Tensor<2,3,double> diffusion_tensor;
        diffusion_tensor[0][0] = 1;
        if constexpr(dim>=2) {
            diffusion_tensor[0][1] = 0.0;
            diffusion_tensor[1][0] = 0.0;
            diffusion_tensor[1][1] = 1.0;
        }
        if constexpr(dim>=3) {
            diffusion_tensor[0][2] = 0.0;
            diffusion_tensor[1][2] = 0.0;
            diffusion_tensor[2][0] = 0.0;
            diffusion_tensor[2][1] = 0.0;
            diffusion_tensor[2][2] = 1.0;
        }
        return diffusion_tensor;
    }

    virtual real objective_function(
        const dealii::Point<dim,real> &pos) const = 0;
};

template <int dim, int nspecies, int nstate, typename real>
class diffusion_u : public diffusion_objective <dim, nspecies, nstate, real>
{
protected:
    // For overloading the virtual functions defined in PhysicsBase
    using Physics::PhysicsBase<dim,nspecies,nstate,real>::dissipative_flux;
    using Physics::PhysicsBase<dim,nspecies,nstate,real>::source_term;
public:
    diffusion_u(
        const Parameters::AllParameters *const parameters_input,
        const bool convection, 
        const bool diffusion): 
            diffusion_objective<dim,nspecies,nstate,real>::diffusion_objective(
                parameters_input,
                convection, 
                diffusion,
                std::make_shared<ManufacturedSolutionU<dim,nspecies,real>>())
    {}

    std::array<real,nstate> source_term (
        const dealii::Point<dim,real> &pos,
        const std::array<real,nstate> &solution,
        const real current_time,
        const dealii::types::global_dof_index cell_index) const override;

    real objective_function(
        const dealii::Point<dim,real> &pos) const override;
};

template <int dim, int nspecies, int nstate, typename real>
class diffusion_v : public diffusion_objective <dim, nspecies, nstate, real>
{
protected:
    // For overloading the virtual functions defined in PhysicsBase
    using Physics::PhysicsBase<dim,nspecies,nstate,real>::dissipative_flux;
    using Physics::PhysicsBase<dim,nspecies,nstate,real>::source_term;
public:
    diffusion_v(
        const Parameters::AllParameters *const parameters_input,
        const bool convection, 
        const bool diffusion): 
            diffusion_objective<dim,nspecies,nstate,real>::diffusion_objective(
                parameters_input,
                convection, 
                diffusion,
                std::make_shared<ManufacturedSolutionV<dim,nspecies,real>>())
    {}

    std::array<real,nstate> source_term (
        const dealii::Point<dim,real> &pos,
        const std::array<real,nstate> &solution,
        const real current_time,
        const dealii::types::global_dof_index cell_index) const override;

    real objective_function(
        const dealii::Point<dim,real> &pos) const override;
};

template <int dim, int nspecies, int nstate, typename real>
class DiffusionFunctional : public Functional<dim, nspecies, nstate, real>
{
    using FadType = Sacado::Fad::DFad<real>; 
    using FadFadType = Sacado::Fad::DFad<FadType>; 
public:
    DiffusionFunctional(
        std::shared_ptr<PHiLiP::DGBase<dim,nspecies,real>> dg_input,
        std::shared_ptr<PHiLiP::Physics::PhysicsBase<dim,nspecies,nstate,FadFadType>> _physics_fad_fad,
        const bool uses_solution_values = true,
        const bool uses_solution_gradient = false)
    : PHiLiP::Functional<dim,nspecies,nstate,real>(dg_input,_physics_fad_fad,uses_solution_values,uses_solution_gradient)
    {}

    template <typename real2>
    real2 evaluate_volume_integrand(
        const PHiLiP::Physics::PhysicsBase<dim,nspecies,nstate,real2> &physics,
        const dealii::Point<dim,real2> &phys_coord,
        const std::array<real2,nstate> &soln_at_q,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &soln_grad_at_q) const;

    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 override
    {
        return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q);
    }
    
    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 override
    {
        return evaluate_volume_integrand<>(physics, phys_coord, soln_at_q, soln_grad_at_q);
    }
};

template <int dim, int nspecies, int nstate>
class DiffusionExactAdjoint : public TestsBase
{
public: 
    explicit DiffusionExactAdjoint(const Parameters::AllParameters *const parameters_input);

    int run_test() const;
};

// for evaluating error slopes
double eval_avg_slope(std::vector<double> error, std::vector<double> grid_size, unsigned int n_grids);

} // Tests namespace
} // PHiLiP namespace

#endif //__DIFFUSION_EXACT_ADJOINT_H__