problem_definition.h

Go to the documentation of this file.

 
#ifndef HYBRIDADRSOLVER_PROBLEM_DEFINITION_H
#define HYBRIDADRSOLVER_PROBLEM_DEFINITION_H
 
#ifndef HYBRIDADRSOLVER_PROBLEM_BASE_H
#define HYBRIDADRSOLVER_PROBLEM_BASE_H
 
#include <cmath>
#include <deal.II/base/function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
#include <set>
#include <string>
 
namespace HybridADRSolver {
using namespace dealii;
using numbers::PI;

template <int dim> class ProblemInterface {
public:
    virtual ~ProblemInterface() = default;

    virtual double diffusion_coefficient(const Point<dim>& p) const = 0;

    virtual Tensor<1, dim> advection_field(const Point<dim>& p) const = 0;

    virtual double reaction_coefficient(const Point<dim>& p) const = 0;

    virtual double source_term(const Point<dim>& p) const = 0;

    [[nodiscard]] virtual std::set<types::boundary_id>
    get_dirichlet_ids() const = 0;

    [[nodiscard]] virtual std::set<types::boundary_id>
    get_neumann_ids() const = 0;

    virtual const Function<dim>&
    get_dirichlet_function(types::boundary_id id) const = 0;

    virtual const Function<dim>&
    get_neumann_function(types::boundary_id id) const = 0;

    [[nodiscard]] virtual bool has_exact_solution() const = 0;

    virtual const Function<dim>& get_exact_solution() const = 0;

    [[nodiscard]] virtual bool is_symmetric() const = 0;

    [[nodiscard]] virtual std::string get_name() const = 0;
};
 
namespace Problems {

template <int dim> class ExactSolution : public Function<dim> {
public:
    double value(const Point<dim>& p, const unsigned int) const override {
        double val = 1.0;
        for (unsigned int d = 0; d < dim; ++d)
            val *= std::sin(PI * p[d]);
        return val;
    }

    Tensor<1, dim> gradient(const Point<dim>& p,
                            const unsigned int) const override {
        Tensor<1, dim> grad;
        for (unsigned int d = 0; d < dim; ++d) {
            grad[d] = PI * std::cos(PI * p[d]);
            for (unsigned int other = 0; other < dim; ++other)
                if (d != other)
                    grad[d] *= std::sin(PI * p[other]);
        }
        return grad;
    }
};

template <int dim> class NeumannFluxRight : public Function<dim> {
public:
    double value(const Point<dim>& p, const unsigned int) const override {
        // Flux = du/dx at x=1.
        // u = sin(pi*x) * sin(pi*y)...
        // du/dx = pi * cos(pi*x) * sin(pi*y)...
        // At x=1, cos(pi*x) = -1.
        double val = -1.0 * PI;
        for (unsigned int d = 1; d < dim; ++d)
            val *= std::sin(PI * p[d]);
        return val;
    }
};

template <int dim> class ADRProblem : public ProblemInterface<dim> {
public:
    ADRProblem() : zero_function(1), neumann_flux(), exact_solution() {}

    double diffusion_coefficient(const Point<dim>&) const override {
        return 1.0;
    }

    double reaction_coefficient(const Point<dim>&) const override {
        return 0.1;
    }

    Tensor<1, dim> advection_field(const Point<dim>& p) const override {
        Tensor<1, dim> beta;
        beta[0] = -p[1];
        beta[1] = p[0];
        if (dim == 3)
            beta[2] = 0.1;
        return beta;
    }

    double source_term(const Point<dim>& p) const override {
        const double u = exact_solution.value(p, 0);
        Tensor<1, dim> grad = exact_solution.gradient(p, 0);
 
        const double mu = diffusion_coefficient(p);
        const double gamma = reaction_coefficient(p);
        Tensor<1, dim> beta = advection_field(p);
 
        // Laplacian of u = product(sin(pi*x_i)) is -d * pi^2 * u
        const double lap = -1.0 * dim * PI * PI * u;
 
        return -mu * lap + beta * grad + gamma * u;
    }

    [[nodiscard]] std::set<types::boundary_id>
    get_dirichlet_ids() const override {
        std::set<types::boundary_id> ids;
        ids.insert(0); // Left
        ids.insert(2); // Bottom
        ids.insert(3); // Top
        if (dim == 3) {
            ids.insert(4);
            ids.insert(5);
        }
        return ids;
    }

    [[nodiscard]] std::set<types::boundary_id>
    get_neumann_ids() const override {
        return {1};
    }

    const Function<dim>&
    get_dirichlet_function(types::boundary_id) const override {
        return zero_function;
    }

    const Function<dim>&
    get_neumann_function(const types::boundary_id id) const override {
        if (id == 1)
            return neumann_flux;
        Assert(false, ExcMessage("Unknown Neumann ID requested"));
        return zero_function;
    }
 
    [[nodiscard]] bool has_exact_solution() const override { return true; }
    const Function<dim>& get_exact_solution() const override {
        return exact_solution;
    }

    [[nodiscard]] bool is_symmetric() const override { return false; }
 
    [[nodiscard]] std::string get_name() const override {
        return "Manufactured ADR (Dirichlet + Neumann)";
    }
 
private:
    Functions::ZeroFunction<dim> zero_function;
    NeumannFluxRight<dim> neumann_flux;
    ExactSolution<dim> exact_solution;
};
 
} // namespace Problems
} // namespace HybridADRSolver
 
#endif
 
#endif // HYBRIDADRSOLVER_PROBLEM_DEFINITION_H