Re: [deal.II] Implementation of a solver for the non-linear diffusion equation using the TimeStepping-Class
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
> void LaplaceOperator::local_apply_cell(
> const MatrixFree & data,
> vector_t & dst,
> const vector_t &src,
> const std::pair & cell_range)
> const
> {
> FEEvaluation > 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
> 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
> void LaplaceOperator::local_apply_cell(
> const MatrixFree & data,
> vector_t & dst,
> const vector_t &src,
> const std::pair & cell_range)
> const
> {
> FEEvaluation > 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
> 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 &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
> inline DEAL_II_ALWAYS_INLINE
> Tensor<1, n_components_to_use, Tensor<1, dim, Number>>
> calculate_gradient_coefficient(
>
[deal.II] Implementation of a solver for the non-linear diffusion equation using the TimeStepping-Class
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
void LaplaceOperator::local_apply_cell(
const MatrixFree & data,
vector_t & dst,
const vector_t &src,
const std::pair & cell_range)
const
{
FEEvaluation
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
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
void LaplaceOperator::local_apply_cell(
const MatrixFree & data,
vector_t & dst,
const vector_t &src,
const std::pair & cell_range)
const
{
FEEvaluation
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
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 &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
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 compon
