VertexCFD_InitialCondition_Gaussian_impl.hpp

#ifndef VERTEXCFD_INITIALCONDITION_GAUSSIAN_IMPL_HPP
#define VERTEXCFD_INITIALCONDITION_GAUSSIAN_IMPL_HPP
 
#include "utils/VertexCFD_Utils_Constants.hpp"
 
#include <Panzer_GlobalIndexer.hpp>
#include <Panzer_HierarchicParallelism.hpp>
#include <Panzer_PureBasis.hpp>
#include <Panzer_Workset_Utilities.hpp>
#include <Teuchos_StandardParameterEntryValidators.hpp>
 
#include <Teuchos_Array.hpp>
 
#include <cmath>
#include <string>
 
namespace VertexCFD
{
namespace InitialCondition
{
//---------------------------------------------------------------------------//
template<class EvalType, class Traits, int NumSpaceDim>
Gaussian<EvalType, Traits, NumSpaceDim>::Gaussian(
    const Teuchos::ParameterList& params, const panzer::PureBasis& basis)
    : _basis_name(basis.name())
    , _scaling(Scaling::unscaled)
{
    const std::string dof_name = params.get<std::string>("Equation Set Name");
    _ic = PHX::MDField<scalar_type, panzer::Cell, panzer::BASIS>(
        dof_name, basis.functional);
    this->addEvaluatedField(_ic);
    this->addUnsharedField(_ic.fieldTag().clone());
    this->setName("Gaussian Initial Condition: " + dof_name);
 
    // Get scaling flag
    if (params.isType<std::string>("Scaling"))
    {
        const auto type_validator = Teuchos::rcp(
            new Teuchos::StringToIntegralParameterEntryValidator<Scaling>(
                Teuchos::tuple<std::string>("scaled", "unscaled"), "unscaled"));
        _scaling = type_validator->getIntegralValue(
            params.get<std::string>("Scaling"));
    }
 
    // Get Gaussian profile parameters
    const auto center = params.get<Teuchos::Array<double>>("Center");
    const auto sigma = params.get<Teuchos::Array<double>>("Sigma");
    _d = params.get<double>("Base");
 
    const double sqrt2pi = std::sqrt(2.0 * Constants::pi);
    for (int dim = 0; dim < basis.dimension(); ++dim)
    {
        double scale;
        if (_scaling == Scaling::scaled)
        {
            scale = params.get<double>("Scale");
        }
        else
        {
            scale = 1.0 / (sqrt2pi * sigma[dim]);
        }
 
        _a[dim] = scale;
        _b[dim] = center[dim];
        _c[dim] = 0.5 / (sigma[dim] * sigma[dim]);
    }
}
 
//---------------------------------------------------------------------------//
template<class EvalType, class Traits, int NumSpaceDim>
void Gaussian<EvalType, Traits, NumSpaceDim>::postRegistrationSetup(
    typename Traits::SetupData sd, PHX::FieldManager<Traits>&)
{
    _basis_index = panzer::getPureBasisIndex(
        _basis_name, (*sd.worksets_)[0], this->wda);
}
 
//---------------------------------------------------------------------------//
template<class EvalType, class Traits, int NumSpaceDim>
void Gaussian<EvalType, Traits, NumSpaceDim>::evaluateFields(
    typename Traits::EvalData workset)
{
    _basis_coords = this->wda(workset).bases[_basis_index]->basis_coordinates;
    auto policy = panzer::HP::inst().teamPolicy<scalar_type, PHX::Device>(
        workset.num_cells);
    Kokkos::parallel_for(this->getName(), policy, *this);
}
 
//---------------------------------------------------------------------------//
template<class EvalType, class Traits, int NumSpaceDim>
KOKKOS_INLINE_FUNCTION void Gaussian<EvalType, Traits, NumSpaceDim>::operator()(
    const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
{
    const int cell = team.league_rank();
    const int num_basis = _ic.extent(1);
 
    Kokkos::parallel_for(
        Kokkos::TeamThreadRange(team, 0, num_basis), [&](const int basis) {
            using Kokkos::exp;
            double result = 1.0;
            for (int dim = 0; dim < num_space_dim; ++dim)
            {
                const auto x = _basis_coords(cell, basis, dim);
                result *= _a[dim]
                          * exp((_b[dim] - x) * (x - _b[dim]) * _c[dim]);
            }
            _ic(cell, basis) = result + _d;
        });
}
 
//---------------------------------------------------------------------------//
 
} // end namespace InitialCondition
} // end namespace VertexCFD
 
#endif // end VERTEXCFD_INITIALCONDITION_GAUSSIAN_IMPL_HPP