Following step-37 I tried to extend it for non-linear equations, using a
jacobian approximation for the "matrix" (which can be written as Jv = (F(u
+ ev) - F(u)) * 1/e, with e a small value (~1e-6), F() the residual of the
function, u the current solution, v the krylov vector obtained from the
GMRES-solver and J the jacobian matrix). Using this method I tried to
calculate the diagonal, which then is used for the preconditioner.
Based on the equation given above and the equation for the linear heat
equation, I get
F(u) = (u - u_0)/dt - (nabla^2u+nabla^2u_0)/2
F(u+ev) = F(u) + ev/dt - e(nabla^2v/2)
Thus Jv can be written as
Jv = (v/dt - nabla^2v/2)
or in weak form
Jv = v\phi/dt + 0.5*\nabla v\nabla\phi
Using that I tried to write the function for creating the diagonal
accordingly:
template <int dim, int fe_degree, typename number>
void JacobianOperator<dim, fe_degree, number>::local_compute_diagonal(
const MatrixFree<dim, number> & data,
vector_t<number> &dst,
const unsigned int &,
const std::pair<unsigned int, unsigned int> &cell_range) const
{
FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(data),
phi_old(data);
AlignedVector<VectorizedArray<number>> diagonal(phi.dofs_per_cell);
for (unsigned int cell = cell_range.first; cell < cell_range.second; ++
cell)
{
phi_old.reinit(cell);
for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < phi.dofs_per_cell; ++j){
phi_old.submit_dof_value(VectorizedArray<number>(), j);
}
phi_old.submit_dof_value(make_vectorized_array<number>(1.), i);
phi_old.evaluate(true, true);
for (unsigned int q = 0; q < phi_old.n_q_points; ++q){
phi_old.submit_gradient(0.5
* (phi_old.get_gradient(q)),
q);
phi_old.submit_value((phi_old.get_value(q) / time_step),
q);
}
phi_old.integrate(true, true);
diagonal[i] = phi_old.get_dof_value(i);
}
for (unsigned int i = 0; i < phi_old.dofs_per_cell; ++i)
phi_old.submit_dof_value(diagonal[i], i);
phi_old.distribute_local_to_global(dst);
}
}
Still, the convergence behaviour is significantly worse than using an
Identity preconditioner. Thus I was wondering if this is related to the
wrong preconditioner, or due to wrong code? I am still trying to understand
how to create the preconditioner for that case most efficiently, thus there
might be some errors here.
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/cd41c84f-9441-4ca9-b48a-d8b9f120456f%40googlegroups.com.
