No, that code is the vmult-call I am using in my LinearOperator (and is
used directly in my solve()-call). The calculate_residual()-function is
similar to the function in step-15, with the difference, that I do not
provide alpha, but the vector, and return the calculated residual value.
Thus, as far as I could understand, I have to do the following in the
matrix-free cell loop:
- Calculate residual of the current solution
- Calculate residual of the current solution + epsilon * src
- Calculate dst, and return it
On the other hand, would it be sufficient to unroll the
residual-calculation, i.e. if the residual is F(u) = nabla^2 u + f, to
write (in the cell-loop): (nabla^2(u + eps*src) + f - nabla^2u -
f)/epsilon? Or would that lead to wrong results?
Am Dienstag, 23. Juli 2019 22:31:06 UTC+2 schrieb Daniel Arndt:
>
> My aim is to implement the jacobian approximation, i.e. (if F(u) = 0 is
>> the function to solve) I wanted to implement
>> dst = (F(current_solution + epsilon * src) - F(current_solution))/epsilon.
>> In my LinearOperator, my code is
>> LinearOperator<LinearAlgebraTrilinos::MPI::Vector> jacobian_operator;
>> jacobian_operator.vmult = [&](LinearAlgebraTrilinos::MPI::Vector &dst,
>> const LinearAlgebraTrilinos::MPI::Vector &src){
>> local_src = src;
>> extended_src = src;
>> local_solution = present_solution;
>>
>> double epsilon = 1e-6;
>> if(local_src.l2_norm() != 0){
>> epsilon = sqrt((1 + local_solution.l2_norm()) * 1e-16) /
>> local_src.l2_norm();
>> };
>> forward_eps_solution = present_solution;
>> backward_eps_solution = present_solution;
>>
>> extended_src *= epsilon;
>> forward_eps_solution += extended_src;
>> backward_eps_solution -= extended_src;
>>
>> compute_residual(backward_eps_solution, cur_residual);
>> compute_residual(forward_eps_solution, eps_residual);
>> eps_residual -= cur_residual;
>> eps_residual /= (2 * epsilon);
>> dst.reinit(locally_owned_dofs,
>> mpi_communicator);
>> dst = eps_residual;
>> };
>>
>> Is the same/similar possible for the matrix-free-approach?
>>
> Sure, you can always write a wrapper around the actual vmult call (which I
> assume is happening in compute_residual()).
>
> Best,
> Daniel
>
>
--
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/5f0807fd-6829-4d60-b49f-41b7e1bb2e27%40googlegroups.com.
