config.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// Copyright (C) 2023 by the Sapphire++ authors
//
// This file is part of Sapphire++.
//
// Sapphire++ is free software: you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License as published by the Free
// Software Foundation, either version 3 of the License, or (at your option) any
// later version.
//
// Sapphire++ is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
// details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with Sapphire++. If not, see <https://www.gnu.org/licenses/>.
//
// -----------------------------------------------------------------------------

 
#ifndef CONFIG_H
#define CONFIG_H
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/function.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/point.h>
 
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/vector.h>
 
#include <cmath>
#include <numbers>
#include <vector>
 
#include "pde-system.h"
#include "sapphirepp-logstream.h"
#include "vfp-flags.h"
 
 
 
namespace sapphirepp
{
  class PhysicalParameters
  {
  public:
    // !!!EDIT HERE!!!
 
 

    PhysicalParameters() = default;
 
 

    void
    declare_parameters(dealii::ParameterHandler &prm)
    {
      dealii::LogStream::Prefix prefix_startup("Startup", saplog);
      dealii::LogStream::Prefix prefix("PhysicalParameters", saplog);
      saplog << "Declaring parameters" << std::endl;
      prm.enter_subsection("Physical parameters");

      // !!!EDIT HERE!!!
 
      prm.leave_subsection();
    }
 
 

    void
    parse_parameters(dealii::ParameterHandler &prm)
    {
      dealii::LogStream::Prefix prefix_startup("Startup", saplog);
      dealii::LogStream::Prefix prefix("PhysicalParameters", saplog);
      saplog << "Parsing parameters" << std::endl;
      prm.enter_subsection("Physical parameters");

      // !!!EDIT HERE!!!
 
      prm.leave_subsection();
    }
  };
 
 
 
  namespace VFP
  {
    // !!!EDIT HERE!!!
    constexpr unsigned int dimension = 2;
 
 

    // !!!EDIT HERE!!!
    constexpr VFPFlags vfp_flags = VFPFlags::time_evolution |          //
                                   VFPFlags::spatial_advection |       //
                                   VFPFlags::momentum |                //
                                   VFPFlags::collision |               //
                                   VFPFlags::source |                  //
                                   VFPFlags::time_independent_fields | //
                                   VFPFlags::time_independent_source;
 
 

    template <unsigned int dim>
    class InitialValueFunction : public dealii::Function<dim>
    {
    public:
      InitialValueFunction(const PhysicalParameters &physical_parameters,
                           const unsigned int        system_size)
        : dealii::Function<dim>(system_size)
        , prm{physical_parameters}
        , lms_indices{PDESystem::create_lms_indices(system_size)}
      {}
 
 

      void
      vector_value([[maybe_unused]] const dealii::Point<dim> &point,
                   dealii::Vector<double>                    &f) const override
      {
        AssertDimension(f.size(), this->n_components);
 
        // NOLINTNEXTLINE(modernize-loop-convert)
        for (unsigned int i = 0; i < f.size(); ++i)
          {
            // !!!EDIT HERE!!!
            f[i] = 0.;
          }
      }
 
 
 
    private:
      const PhysicalParameters prm;
      const std::vector<std::array<unsigned int, 3>> lms_indices;
    };
 
 

    template <unsigned int dim>
    class BoundaryValueFunction : public dealii::Function<dim>
    {
    public:
 
      BoundaryValueFunction(const PhysicalParameters &physical_parameters,
                            const unsigned int        system_size)
        : dealii::Function<dim>(system_size)
        , prm{physical_parameters}
        , lms_indices{PDESystem::create_lms_indices(system_size)}
      {}
 
 

      void
      bc_vector_value_list(
        [[maybe_unused]] const std::vector<dealii::Point<dim>> &points,
        [[maybe_unused]] const unsigned int                     boundary_id,
        [[maybe_unused]] std::vector<dealii::Vector<double>>   &bc_values) const
      {
        AssertDimension(points.size(), bc_values.size());
        AssertDimension(bc_values[0].size(), this->n_components);
 
        for (unsigned int q_index = 0; q_index < points.size(); ++q_index)
          {
            // !!!EDIT HERE!!!
            if (boundary_id == 0)
              {
                // lower x
              }
            else if (boundary_id == 1)
              {
                // upper x
              }
            else if (boundary_id == 2)
              {
                // lower p
              }
            else if (boundary_id == 3)
              {
                // upper p
              }
          }
      }
 
 
 
    private:
      const PhysicalParameters prm;
      const std::vector<std::array<unsigned int, 3>> lms_indices;
    };
 
 

    template <unsigned int dim>
    class ScatteringFrequency : public dealii::Function<dim>
    {
    public:
      ScatteringFrequency(const PhysicalParameters &physical_parameters)
        : dealii::Function<dim>(1)
        , prm{physical_parameters}
      {}
 
 

      void
      value_list(
        const std::vector<dealii::Point<dim>> &points,
        std::vector<double>                   &scattering_frequencies,
        [[maybe_unused]] const unsigned int    component = 0) const override
      {
        AssertDimension(scattering_frequencies.size(), points.size());
 
        for (unsigned int q_index = 0; q_index < points.size(); ++q_index)
          {
            // !!!EDIT HERE!!!
            scattering_frequencies[q_index] = 100.;
          }
      }
 
 
 
    private:
      const PhysicalParameters prm;
    };
 
 

    template <unsigned int dim>
    class Source : public dealii::Function<dim>
    {
    public:
      Source(const PhysicalParameters &physical_parameters,
             unsigned int              system_size)
        : dealii::Function<dim>(system_size)
        , prm{physical_parameters}
        , lms_indices{PDESystem::create_lms_indices(system_size)}
      {}
 

      void
      vector_value([[maybe_unused]] const dealii::Point<dim> &point,
                   dealii::Vector<double> &source_values) const override
      {
        AssertDimension(source_values.size(), this->n_components);
 
        // NOLINTNEXTLINE(modernize-loop-convert)
        for (unsigned int i = 0; i < source_values.size(); ++i)
          {
            // !!!EDIT HERE!!!
            // Define parameters
            const double p0    = 1.;
            const double Q     = 0.1;
            const double sig_p = 0.1;
            const double sig_x = 0.001;
 
            // Get x and p value
            const double x = point[0];
            const double p = std::exp(point[1]); // point[1] = log(p)
 
            // Convert i -> l, m, s
            const unsigned int l = lms_indices[i][0];
            const unsigned int m = lms_indices[i][1];
            const unsigned int s = lms_indices[i][2];
 
            // isotropic part
            if (l == 0 && m == 0 && s == 0)
              {
                // shifted Gaussian in x and p
                // s_000 = sqrt(4 pi) * s
                source_values[i] =
                  Q / (std::sqrt(std::numbers::pi) * sig_p * sig_x) *
                  std::exp(-(p - p0) * (p - p0) / (2. * sig_p * sig_p)) *
                  std::exp(-x * x / (2. * sig_x * sig_x));
              }
            // vanishing anisotropic part
            else
              {
                source_values[i] = 0.;
              }
          }
      }
 
 
 
    private:
      const PhysicalParameters prm;
      const std::vector<std::array<unsigned int, 3>> lms_indices;
    };
 
 

    template <unsigned int dim>
    class MagneticField : public dealii::Function<dim>
    {
    public:
      MagneticField(const PhysicalParameters &physical_parameters)
        : dealii::Function<dim>(3)
        , prm{physical_parameters}
      {}
 

      void
      vector_value([[maybe_unused]] const dealii::Point<dim> &point,
                   dealii::Vector<double> &magnetic_field) const override
      {
        AssertDimension(magnetic_field.size(), this->n_components);

        // !!!EDIT HERE!!!
        magnetic_field[0] = 0.; // B_x
        magnetic_field[1] = 0.; // B_y
        magnetic_field[2] = 0.; // B_z
      }
 
 
 
    private:
      const PhysicalParameters prm;
    };
 
 

    template <unsigned int dim>
    class BackgroundVelocityField : public dealii::Function<dim>
    {
    public:
      BackgroundVelocityField(const PhysicalParameters &physical_parameters)
        : dealii::Function<dim>(3)
        , prm{physical_parameters}
      {}
 
 

      void
      vector_value([[maybe_unused]] const dealii::Point<dim> &point,
                   dealii::Vector<double> &velocity) const override
      {
        AssertDimension(velocity.size(), this->n_components);

        // !!!EDIT HERE!!!
        // u(x) = u_sh/2r * ((1-r)*tanh(x/d_sh) + (1+r))
        const double r    = 4.;
        const double d_sh = 0.001;
        const double u_sh = 0.1;
 
        // u_x
        velocity[0] =
          u_sh / (2 * r) * ((1 - r) * std::tanh(point[0] / d_sh) + (1 + r));
        velocity[1] = 0.; // u_y
        velocity[2] = 0.; // u_z
      }
 

      void
      divergence_list(
        [[maybe_unused]] const std::vector<dealii::Point<dim>> &points,
        std::vector<double> &divergence) const
      {
        AssertDimension(divergence.size(), points.size());
 
        for (unsigned int q_index = 0; q_index < points.size(); ++q_index)
          {
            // !!!EDIT HERE!!!
            // u(x) = u_sh/2r * ((1-r)*tanh(x/d_sh) + (1+r))
            // => d/dx u(x) = u_sh/2r 1/d_sh (1-r) (1-tanh(x/d_sh)^2)
            const double r    = 4.;
            const double d_sh = 0.001;
            const double u_sh = 0.1;
 
            const double x = points[q_index][0];
 
            // div u
            divergence[q_index] =
              u_sh / (2 * r) * (1 - r) / d_sh *
              (1 - std::tanh(x / d_sh) * std::tanh(x / d_sh));
          }
      }
 
 

      void
      material_derivative_list(
        const std::vector<dealii::Point<dim>> &points,
        std::vector<dealii::Vector<double>>   &material_derivatives) const
      {
        AssertDimension(material_derivatives.size(), points.size());
        AssertDimension(material_derivatives[0].size(), this->n_components);
 
        for (unsigned int q_index = 0; q_index < points.size(); ++q_index)
          {
            // !!!EDIT HERE!!!
            // u(x) = u_sh/2r * ((1-r)*tanh(x/d_sh) + (1+r))
            // => D/Dt u(x) = d/dt u(x) + u d/dx u(x)
            //     = (u_sh/2r)^2 * ((1-r)*tanh(x/d_sh) + (1+r)) *
            //        1/d_sh (1-r) (1-tanh(x/d_sh)^2)
            const double r    = 4.;
            const double d_sh = 0.001;
            const double u_sh = 0.1;
 
            const double x = points[q_index][0];
 
            // D/Dt u_x
            material_derivatives[q_index][0] =
              u_sh * u_sh / (4 * r * r) / d_sh * (1 - r) *
              ((1 - r) * std::tanh(x / d_sh) + (1 + r)) *
              (1 - std::tanh(x / d_sh) * std::tanh(x / d_sh));
            material_derivatives[q_index][1] = 0.; // D/Dt u_y
            material_derivatives[q_index][2] = 0.; // D/Dt u_z
          }
      }
 
 

      void
      jacobian_list(const std::vector<dealii::Point<dim>>   &points,
                    std::vector<dealii::FullMatrix<double>> &jacobians) const
      {
        AssertDimension(jacobians.size(), points.size());
        AssertDimension(jacobians[0].m(), this->n_components);
        AssertDimension(jacobians[0].n(), this->n_components);
 
        for (unsigned int q_index = 0; q_index < points.size(); ++q_index)
          {
            // !!!EDIT HERE!!!
            //  u(x) = u_sh/2r * ((1-r)*tanh(x/d_sh) + (1+r))
            //  => u_00 = du/dx = u_sh/2r 1/d_sh (1-r) (1-tanh(x/d_sh)^2)
            const double r    = 4.;
            const double d_sh = 0.001;
            const double u_sh = 0.1;
 
            const double x = points[q_index][0];
 
            // \partial u_x / \partial x
            jacobians[q_index][0][0] =
              u_sh / (2 * r) * (1 - r) / d_sh *
              (1 - std::tanh(x / d_sh) * std::tanh(x / d_sh));
            jacobians[q_index][0][1] = 0.; // \partial u_x / \partial y
            jacobians[q_index][0][2] = 0.; // \partial u_x / \partial z
 
            jacobians[q_index][1][0] = 0.; // \partial u_y / \partial x
            jacobians[q_index][1][1] = 0.; // \partial u_y / \partial y
            jacobians[q_index][1][2] = 0.; // \partial u_y / \partial z
 
            jacobians[q_index][2][0] = 0.; // \partial u_z / \partial x
            jacobians[q_index][2][1] = 0.; // \partial u_z / \partial y
            jacobians[q_index][2][2] = 0.; // \partial u_z / \partial z
          }
      }
 
 
 
    private:
      const PhysicalParameters prm;
    };
  } // namespace VFP
} // namespace sapphirepp
#endif