[petsc-users] Approximation of matrix vector product by finite difference method

Fri Sep 7 08:16:45 CDT 2018

On Fri, Sep 7, 2018 at 6:39 AM Yingjie Wu <yjwu16 at gmail.com> wrote:

> Thank you very much for your reply.
>
> I'm a bit confused if I use what you recommend:  - snes_fd -pc_type lu
> means that I explicitly construct the Jacobian matrix using the finite
> difference method, construct the precondition matrix using the completely
> LU decomposition, and solve the step size \ delta x with the GMRES method
> (default).
>

You should do this first in order to check that your code reproduces the
right physics.

Then you should do -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type
greedy for larger problems.

Finally, if you want to go even larger, then you use -snes_mf and a user
preconditioner. However, there is not
much benefit of MF with a low order method because you will have to form
some preconditioner.

  Thanks,

     Matt

> In fact, what I want to use is to approximate the vector product of a
> matrix with finite difference, so that the explicit construction of
> Jacobian matrices can be avoided. If so, should I use MatrixFreeMethod? How
> should I set it up? If I want to set up precondition, what do I need to
> add?
>
> In addition, I want to output variables in each nolinear step. What should
> I add code to make SNES step by step?
> There may be many problems, but they bother me very much. I am looking
> forward to your reply.
>
> Thanks,
> Yingjie
>
>
> Matthew Knepley <knepley at gmail.com> 于2018年9月6日周四 下午10:34写道：
>
>> On Thu, Sep 6, 2018 at 4:47 AM Yingjie Wu <yjwu16 at gmail.com> wrote:
>>
>>> Dear Petsc developer:
>>> Hi,
>>> Thank you for your previous help.
>>> I recently modeled on PETSc's SNES example and wrote a computer program
>>> myself. This program is mainly for solving nonlinear equations of thermal
>>> hydraulics.
>>>
>>> ∇·（λ∇T） - ∇_y（ρ*Cp*u） - T_source = 0
>>>
>>> w*ρ*u = ρg - ∇_y（P）
>>>
>>>      ∇·（ 1/w * ∇P ） =  - ∇（ ρg / w ）
>>>
>>> Where P, T and u are variables, the distribution represents pressure,
>>> temperature and velocity. The rest are nonlinear physical parameters and
>>> constants.
>>>
>>> Because the program is very preliminary, so I use - snes_mf so that I
>>> can save the part of writing to calculate the Jacobian matrix.
>>>
>>> After compiling and passing, I found that the residual function had not
>>> dropped to a small enough level, but the program stopped, as follows:
>>>
>>> Setting Up: -snes_mf  -snes_monitor  -draw_pause  10  -snes_view
>>>
>>>
>> First, do not use -snes_mf. It is not for testing, but for sophisticated
>> use. The first option you
>> might try is
>>
>>    -snes_fd -pc_type lu
>>
>> That uses a full Jacobian and LU factorization for a direct solve. Always
>> run the solve using
>>
>>   -snes_view -snes_converged_reason -snes_monitor -ksp_converged_reason
>> -ksp_monitor_true_residual
>>
>> When that gets too expensive, you can try
>>
>>   -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy
>>
>> but that requires you to preallocate the Jacobian matrix correctly.
>>
>>   Thanks,
>>
>>     Matt
>>
>>> 0 SNES Function norm 3.724996516631e+09
>>>
>>>   1 SNES Function norm 2.194322909557e+09
>>>
>>>   2 SNES Function norm 1.352051559826e+09
>>>
>>>   3 SNES Function norm 1.522311916217e+08
>>>
>>> SNES Object: 1 MPI processes
>>>
>>>   type: newtonls
>>>
>>>   maximum iterations=50, maximum function evaluations=10000
>>>
>>>   tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
>>>
>>>   total number of linear solver iterations=1298
>>>
>>>   total number of function evaluations=11679
>>>
>>>   norm schedule ALWAYS
>>>
>>>   SNESLineSearch Object: 1 MPI processes
>>>
>>>     type: bt
>>>
>>>       interpolation: cubic
>>>
>>>       alpha=1.000000e-04
>>>
>>>     maxstep=1.000000e+08, minlambda=1.000000e-12
>>>
>>>     tolerances: relative=1.000000e-08, absolute=1.000000e-15,
>>> lambda=1.000000e-08
>>>
>>>     maximum iterations=40
>>>
>>>   KSP Object: 1 MPI processes
>>>
>>>     type: gmres
>>>
>>>       restart=30, using Classical (unmodified) Gram-Schmidt
>>> Orthogonalization with no iterative refinement
>>>
>>>       happy breakdown tolerance 1e-30
>>>
>>>     maximum iterations=10000, initial guess is zero
>>>
>>>     tolerances:  relative=1e-05, absolute=1e-50, divergence=10000.
>>>
>>>     left preconditioning
>>>
>>>     using PRECONDITIONED norm type for convergence test
>>>
>>>   PC Object: 1 MPI processes
>>>
>>>     type: none
>>>
>>>     linear system matrix = precond matrix:
>>>
>>>     Mat Object: 1 MPI processes
>>>
>>>       type: mffd
>>>
>>>       rows=300, cols=300
>>>
>>>         Matrix-free approximation:
>>>
>>>           err=1.49012e-08 (relative error in function evaluation)
>>>
>>>           Using wp compute h routine
>>>
>>>               Does not compute normU
>>>
>>> I would like to know why the residual function can not continue to
>>> decline, and why the program will stop before convergence.
>>> I do not know much about the convergence criteria and convergence rules
>>> of PETSc for solving nonlinear equations. I hope I can get your help.
>>> I'm looking forward to your reply~
>>>
>>> Thanks,
>>> Yingjie
>>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> <http://www.cse.buffalo.edu/~knepley/>
>>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180907/a25d2cab/attachment-0001.html>