Found my mistake, I was iterating over the wrong things.
Am Sonntag, 19. Januar 2020 14:54:16 UTC+1 schrieb 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!
>
> Am Samstag, 18. Januar 2020 12:11:43 UTC+1 schrieb 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!
>>
>
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/b1dc11b3-01b7-40a8-953f-5d694469c1b2%40googlegroups.com.
