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,


On Sat, 18 Jan 2020, 14:41 'Maxi Miller' via deal.II User Group, <
dealii@googlegroups.com> wrote:

> 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!
>
> --
> 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+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/615d71a1-fb27-4769-8766-5f17872a9053%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/615d71a1-fb27-4769-8766-5f17872a9053%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/CAKGOWiN1%2B6FH7yH2WfyssYAfSaSnL5uQhUgoUUO%2Boq6yKMyEKg%40mail.gmail.com.

Reply via email to