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/615d71a1-fb27-4769-8766-5f17872a9053%40googlegroups.com.
