Implementation of a solver for the non-linear diffusion equation using the TimeStepping-Class

[Maxi Miller]

I tried to implement a solver for the non-linear diffusion equation (\partial_t u = grad(u(grad u)) - f) using the TimeStepping-Class, the EmbeddedExplicitRungeKutta-Method and (for assembly) the matrix-free approach. For initial tests I used the linear heat equation with the solution u = exp(-pi * pi * t) * sin(pi * x) * sin(pi * y) and the assembly routine   

template <int dim, int degree, int n_points_1d>
    void LaplaceOperator<dim, degree, n_points_1d>::local_apply_cell(
            const MatrixFree<dim, Number> &                   data,
            vector_t &      dst,
            const vector_t &src,
            const std::pair<unsigned int, unsigned int> &     cell_range) const
    {
        FEEvaluation<dim, degree, n_points_1d, n_components_to_use, Number> phi(data);

        for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
        {
            phi.reinit(cell);
            phi.read_dof_values_plain(src);
            phi.evaluate(false, true);
            for (unsigned int q = 0; q < phi.n_q_points; ++q)
            {
                auto gradient_coefficient = calculate_gradient_coefficient(phi.get_gradient(q));
                phi.submit_gradient(gradient_coefficient, q);
            }
            phi.integrate_scatter(false, true, dst);
        }
    }

and

    template <int dim, typename Number>
    inline DEAL_II_ALWAYS_INLINE
    Tensor<1, n_components_to_use, Tensor<1, dim, Number>> calculate_gradient_coefficient(
        #if defined(USE_NONLINEAR) || defined(USE_ADVECTION)
            const Tensor<1, n_components_to_use, Number> &input_value,
        #endif
            const Tensor<1, n_components_to_use, Tensor<1, dim, Number>> &input_gradient){
        Tensor<1, n_components_to_use, Tensor<1, dim, Number>> ret_val;
        for(size_t component = 0; component < n_components_to_use; ++component){
            for(size_t d = 0; d < dim; ++d){
                ret_val[component][d] = -1. / (2 * M_PI * M_PI) * input_gradient[component][d];
            }
        }
        return ret_val;
    }

This approach works, and delivers correct results. Now I wanted to test the same approach for the non-linear diffusion equation with f = -exp(-2 * pi^2 * t) * 0.5 * pi^2 * (-cos(2 * pi * (x - y)) - cos(2 * pi * (x + y)) + cos(2 * pi * x) + cos(2 * pi * y)), which should be the solution to grad(u (grad u)) with u = exp(-pi^2*t) * sin(pi * x) * sin(pi * y). Thus, I changed the routines to

    template <int dim, int degree, int n_points_1d>
    void LaplaceOperator<dim, degree, n_points_1d>::local_apply_cell(
            const MatrixFree<dim, Number> &                   data,
            vector_t &      dst,
            const vector_t &src,
            const std::pair<unsigned int, unsigned int> &     cell_range) const
    {
        FEEvaluation<dim, degree, n_points_1d, n_components_to_use, Number> phi(data);

        for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
        {
            phi.reinit(cell);
            phi.read_dof_values_plain(src);
            phi.evaluate(true, true);
            for(size_t q = 0; q < phi.n_q_points; ++q){
                auto value = phi.get_value(q);
                auto gradient = phi.get_gradient(q);
                phi.submit_value(calculate_value_coefficient(value,
                                                             phi.quadrature_point(q),
                                                             local_time), q);
                phi.submit_gradient(calculate_gradient_coefficient(value,
                                                                   gradient), q);
            }
            phi.integrate_scatter(true, true, dst);
        }
    }

and

    template <int dim, typename Number>
    inline DEAL_II_ALWAYS_INLINE
    Tensor<1, n_components_to_use, Number> calculate_value_coefficient(const Tensor<1, n_components_to_use, Number> &input_value,
                                                                       const Point<dim, Number> &point,
                                                                       const double &time){
        Tensor<1, n_components_to_use, Number> ret_val;
        //(void) input_value;
        (void) input_value;
        for(size_t component = 0; component < n_components_to_use; ++component){
            const double x = point[component][0];
            const double y = point[component][1];
            ret_val[component] = (- exp(-2 * M_PI * M_PI * time)
                                  * 0.5 * M_PI * M_PI * (-cos(2 * M_PI * (x - y))
                                                         - cos(2 * M_PI * (x + y))
                                                         + cos(2 * M_PI * x)
                                                         + cos(2 * M_PI * y)))
                    ;
        }
        return ret_val;
    }

    template <int dim, typename Number>
    inline DEAL_II_ALWAYS_INLINE
    Tensor<1, n_components_to_use, Tensor<1, dim, Number>> calculate_gradient_coefficient(
        #if defined(USE_NONLINEAR) || defined(USE_ADVECTION)
            const Tensor<1, n_components_to_use, Number> &input_value,
        #endif
            const Tensor<1, n_components_to_use, Tensor<1, dim, Number>> &input_gradient){
        Tensor<1, n_components_to_use, Tensor<1, dim, Number>> ret_val;
        for(size_t component = 0; component < n_components_to_use; ++component){
            for(size_t d = 0; d < dim; ++d){
                ret_val[component][d] = -1. * input_value[component] * input_gradient[component][d];
            }
        }
        return ret_val;
    }

Unfortunately, now the result is not correct anymore. The initial sin-shape is spreading out in both x- and y-direction, leading to wrong results.Therefore I was wondering if I just implemented the functions in a wrong way, or if there could be something else wrong?

Thanks!

[masod sadipour]

Dear Maxi,

I am not an expert in deal.ii. I have a project which is very similar to your project. I was wondering if it is possible to contact you.

Best,

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/615d71a1-fb27-4769-8766-5f17872a9053%40googlegroups.com.

[Maxi Miller]

I added a small MWE-code to the question, to make debugging easier. Switching between the linear heat equation and the non-linear heat equation can be done using the constexpr-variable use_nonlinear in line 79.

Furthermore, I changed the solution from exp(-pi^2*t) to exp(-2 * pi^2 * t), and changed the value of f accordingly.

Thanks!

[Maxi Miller]

Found my mistake, I was iterating over the wrong things.