[petsc-users] Question about using MatSNESMFWPSetComputeNormU
Choi Kyungjun
kyungjun.choi92 at gmail.com
Thu Oct 6 09:23:48 CDT 2016
Dear Matt.
Thank you very much for your help.
I'm currently working on PETSc library into 2-D compressible Euler / NS
equation solver, especially for convergence of steady state problem.
I adjusted my flow code as you told me, using snes_mf command line option,
but I have a question about SNESsetFunction, especially function evaluation
routine.
My command line option goes like this, as you told me.
*-snes_mf -pc_type none -snes_view -snes_monitor -ksp_monitor
-snes_converged_reason -ksp_converged_reason*
I remember that if I use snes_mf option, the matrix-free method is applied
with computing Jacobian like below.
[image: 본문 이미지 1](captured from Petsc Manual p.113)
But I computed Jacobian with function evaluation routine, with
SNESSetFunction(snes, r, FormPetscResidual, userctx, ier).
I referred to my reference code which computes Jacobian like below.
F'(u) a = -F(u) - (volume)/dt *a
This is just reverse calculation of equation, not matrix-free form. This is
done at the function evaluation routine (FormPetscResidual).
*I want to ask how I can use the REAL matrix-free form.*
I'll attach my flow code and computation log below.
Thank you so much every time for your sincere help.
Kyungjun.
*================================================================================*
*(This is the flow code, and vectors are already created.)*
*call SNESCreate(PETSC_COMM_WORLD, Mixt%snes, ier)*
*call SNESSetFunction(Mixt%snes, Mixt%r, FormPetscResidual, Collect, ier)*
*call SNESSetFromOptions(Mixt%snes, ier)*
*================================================================================*
*(This is function evaluation routine)*
*subroutine FormPetscResidual(snes, x, f, Collect, ier)*
type(t_Collect), intent(inout) :: Collect
SNES :: snes
Vec :: x, f
integer :: ier, counter, iCell, iVar, temp
integer :: ndof
real(8), allocatable :: CVar(:,:)
real(8), allocatable :: PVar(:,:)
PetscScalar, pointer :: xx_v(:)
PetscScalar, pointer :: ff_v(:)
! Set degree of freedom of this system.
ndof = Collect%pMixt%nCVar * Collect%pGrid%nCell
! Backup the original values for cv to local array CVar
allocate( CVar(0:Collect%pMixt%nCVar-1, Collect%pGrid%nCell) )
allocate( PVar(0:Collect%pMixt%nPVar-1, Collect%pGrid%nCell) )
allocate( xx_v(1:ndof) )
allocate( ff_v(1:ndof) )
xx_v(:) = 0d0
ff_v(:) = 0d0
! Backup the original values for cv and pv
do iCell = 1, Collect%pGrid%nCell
do iVar = 0, Collect%pMixt%nCVar-1
CVar(iVar,iCell) = Collect%pMixt%cv(iVar,iCell)
PVar(iVar,iCell) = Collect%pMixt%pv(iVar,iCell)
end do
end do
! Copy the input argument vector x to array value xx_v
call VecGetArrayReadF90(x, xx_v, ier)
call VecGetArrayF90(f, ff_v, ier)
! Compute copy the given vector into Mixt%cv and check for validity
counter = 0
do iCell = 1, Collect%pGrid%nCell
do iVar = 0, Collect%pMixt%nCVar-1
counter = counter + 1
Collect%pMixt%cv(iVar,iCell) = xx_v(counter)
end do
end do
! Update primitive variables with input x vector to compute residual
call PostProcessing(Collect%pMixt,Collect%pGrid,Collect%pConf)
! Compute the residual
call ComputeResidual(Collect%pMixt,Collect%pGrid,Collect%pConf) -->
where update residual of cell
! Copy the residual array into the PETSc vector
counter = 0
do iCell = 1, Collect%pGrid%nCell
do iVar = 0, Collect%pMixt%nCVar-1
counter = counter + 1
* ff_v(counter) = Collect%pMixt%Residual(iVar,iCell) +
Collect%pGrid%vol(iCell)/Collect%pMixt%TimeStep(iCell)*(
Collect%pMixt%cv(iVar,iCell) - CVar(iVar,iCell) )*
end do
end do
! Restore conservative variables
do iCell = 1, Collect%pGrid%nCell
do iVar = 0, Collect%pMixt%nCVar-1
Collect%pMixt%cv(iVar,iCell) = CVar(iVar,iCell)
Collect%pMixt%pv(iVar,iCell) = PVar(iVar,iCell)
end do
end do
call VecRestoreArrayReadF90(x, xx_v, ier)
call VecRestoreArrayF90(f, ff_v, ier)
deallocate(CVar)
deallocate(PVar)
*end subroutine*
*================================================================================*
*Computation log*
[image: 본문 이미지 2]
2016-08-19 21:14 GMT+09:00 Barry Smith <bsmith at mcs.anl.gov>:
>
> It looks like the SNESView() you have below was called before you ever
> did a solve, hence it prints the message "information may be incomplete".
> Note also zero function evaluations have been done in the SNESSolve, if the
> solve had been called it should be great than 0.
>
> SNES Object: 1 MPI processes
> type: newtonls
> SNES has not been set up so information may be incomplete
>
> This is also why it prints
>
> The compute h routine has not yet been set
>
> The information about the h routine won't be printed until after an actual
> solve is done and the "compute h" function is set.
>
> Barry
>
> Note you can call MatMFFDSetType() to control the "compute h" function
> that is used.
>
>
>
> > On Aug 19, 2016, at 12:04 AM, 최경준 <kyungjun.choi92 at gmail.com> wrote:
> >
> > Dear Barry and Matt.
> >
> > Thank you very much for helping me up all night. (in my time)
> >
> > And sorry for not asking with sufficient source code condition or my
> circumstances. (also with poor English.)
> >
> >
> > I just want to make sure that the options of my code is well applied.
> >
> > I'm trying to use GMRES with matrix-free method. I'd like to solve 2-D
> euler equation without preconditioning matrix, for now.
> >
> >
> > 1) I'm still curious whether my snes context is using MF jacobian. (
> just like -snes_mf command line option)
> >
> > 2) And mind if I ask you that whether I applied petsc functions properly?
> >
> > I'll check out ex5 for applying command line options.
> >
> >
> > I'll attach my petsc flow code and option log by SNESView() below.
> > ------------------------------------------------------------
> --------------------------------------------------------
> > - petsc flow code
> > ------------------------------------------------------------
> --------------------------------------------------------
> >
> > ndof = Mixt%nCVar * Grid%nCell
> >
> > call VecCreateMPIWIthArray(PETSC_COMM_WORLD, Mixt%nCVar, ndof,
> PETSC_DECIDE, Mixt%cv, Mixt%x, ier)
> > call VecDuplicate(Mixt%x, Mixt%r, ier)
> > call VecSet(Mixt%r, zero, ier)
> >
> > call SNESCreate(PETSC_COMM_WORLD, Mixt%snes, ier)
> > call SNESSetFunction(Mixt%snes, Mixt%r, FormPetscResidual, Collect, ier)
> > call MatCreateSNESMF(Mixt%snes, Mixt%A, ier)
> >
> > call SNESSetJacobian(Mixt%snes, Mixt%A, Mixt%A, MatMFFDComputeJacobian,
> Collect, ier)
> > call SNESSetFromOptions(Mixt%snes, ier)
> >
> > call SNESGetKSP(Mixt%snes, ksp, ier)
> > call KSPSetType(ksp, KSPGMRES, ier)
> > call KSPGetPC(ksp, pc, ier)
> > call PCSetType(pc, PCNONE, ier)
> > call KSPSetInitialGuessNonzero(ksp, PETSC_TRUE, ier)
> > call KSPGMRESSetRestart(ksp, 30, ier)
> > call KSPGMRESSetPreAllocation(ksp, ier)
> >
> >
> > call SNESSetFunction(Mixt%snes, Mixt%r, FormPetscResidual, Collect, ier)
> > call SNESSetJacobian(Mixt%snes, Mixt%A, Mixt%A, MatMFFDComputeJacobian,
> Collect, ier)
> >
> > call SNESSolve(Mixt%snes, PETSC_NULL_OBJECT, Mixt%x, ier)
> >
> > stop ( for temporary )
> >
> >
> > ------------------------------------------------------------
> --------------------------------------------------------
> > subroutine FormPetscResidual(snes, x, f, Collect, ier)
> > type(t_Collect), intent(inout) :: Collect
> >
> > SNES :: snes
> > Vec :: x, f
> > integer :: ier, counter, iCell, iVar, temp
> > integer :: ndof
> > real(8), allocatable :: CVar(:,:)
> > real(8), allocatable :: PVar(:,:)
> > PetscScalar, pointer :: xx_v(:)
> > PetscScalar, pointer :: ff_v(:)
> >
> > ! Set degree of freedom of this system.
> > ndof = Collect%pMixt%nCVar * Collect%pGrid%nCell
> >
> > ! Backup the original values for cv to local array CVar
> > allocate( CVar(0:Collect%pMixt%nCVar-1, Collect%pGrid%nCell) )
> > allocate( PVar(0:Collect%pMixt%nPVar-1, Collect%pGrid%nCell) )
> > allocate( xx_v(1:ndof) )
> > allocate( ff_v(1:ndof) )
> > xx_v(:) = 0d0
> > ff_v(:) = 0d0
> >
> > ! Backup the original values for cv and pv
> > do iCell = 1, Collect%pGrid%nCell
> > do iVar = 0, Collect%pMixt%nCVar-1
> > CVar(iVar,iCell) = Collect%pMixt%cv(iVar,iCell)
> > PVar(iVar,iCell) = Collect%pMixt%pv(iVar,iCell)
> > end do
> > end do
> >
> > ! Copy the input argument vector x to array value xx_v
> > call VecGetArrayReadF90(x, xx_v, ier)
> > call VecGetArrayF90(f, ff_v, ier)
> >
> > ! Compute copy the given vector into Mixt%cv and check for validity
> > counter = 0
> > do iCell = 1, Collect%pGrid%nCell
> > do iVar = 0, Collect%pMixt%nCVar-1
> > counter = counter + 1
> > Collect%pMixt%cv(iVar,iCell) = xx_v(counter)
> > end do
> > end do
> >
> > ! Update primitive variables with input x vector to compute residual
> > call PostProcessing(Collect%pMixt,Collect%pGrid,Collect%pConf)
> >
> >
> > ! Compute the residual
> > call ComputeResidual(Collect%pMixt,Collect%pGrid,Collect%pConf) -->
> where update residual of cell
> >
> > ! Copy the residual array into the PETSc vector
> > counter = 0
> > do iCell = 1, Collect%pGrid%nCell
> > do iVar = 0, Collect%pMixt%nCVar-1
> > counter = counter + 1
> >
> > ff_v(counter) = Collect%pMixt%Residual(iVar,iCell) +
> Collect%pGrid%vol(iCell)/Collect%pMixt%TimeStep(iCell)*(
> Collect%pMixt%cv(iVar,iCell) - CVar(iVar,iCell) )
> > end do
> > end do
> >
> > ! Restore conservative variables
> > do iCell = 1, Collect%pGrid%nCell
> > do iVar = 0, Collect%pMixt%nCVar-1
> > Collect%pMixt%cv(iVar,iCell) = CVar(iVar,iCell)
> > Collect%pMixt%pv(iVar,iCell) = PVar(iVar,iCell)
> > end do
> > end do
> >
> > call VecRestoreArrayReadF90(x, xx_v, ier)
> > call VecRestoreArrayF90(f, ff_v, ier)
> >
> > deallocate(CVar)
> > deallocate(PVar)
> > ------------------------------------------------------------
> --------------------------------------------------------
> >
> >
> > ------------------------------------------------------------
> --------------------------------------------------------
> > - option log
> > ------------------------------------------------------------
> --------------------------------------------------------
> > SNES Object: 1 MPI processes
> > type: newtonls
> > SNES has not been set up so information may be incomplete
> > maximum iterations=1, maximum function evaluations=10000
> > tolerances: relative=1e-08, absolute=1e-32, solution=1e-08
> > total number of linear solver iterations=0
> > total number of function evaluations=0
> > 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
> > GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
> Orthogonalization with no iterative refinement
> > GMRES: happy breakdown tolerance 1e-30
> > maximum iterations=10000
> > tolerances: relative=0.001, absolute=1e-50, divergence=10000.
> > left preconditioning
> > using nonzero initial guess
> > using DEFAULT norm type for convergence test
> > PC Object: 1 MPI processes
> > type: none
> > PC has not been set up so information may be incomplete
> > linear system matrix = precond matrix:
> > Mat Object: 1 MPI processes
> > type: mffd
> > rows=11616, cols=11616
> > Matrix-free approximation:
> > err=1.49012e-08 (relative error in function evaluation)
> > The compute h routine has not yet been set
> >
> >
> > Sincerely,
> >
> > Kyungjun
> >
> >
> > 2016-08-19 13:00 GMT+09:00 Barry Smith <bsmith at mcs.anl.gov>:
> >
> > > On Aug 18, 2016, at 10:28 PM, 최경준 <kyungjun.choi92 at gmail.com> wrote:
> > >
> > > Dear Matt.
> > >
> > > I didn't use the command line options because it looked not working.
> > >
> > > I called SNESSetFromOptions(snes, ier) in my source code,
> > >
> > > but options like -snes_mf or -snes_monitor doesn't look working.
> >
> > "doesn't work" is not useful to help us figure out what has gone
> wrong. You need to show us EXACTLY what you did by sending the code you
> compiled and the command line options you ran and all the output include
> full error messages. Without the information we simply do not have enough
> information to even begin to guess why it "doesn't work".
> >
> > Barry
> >
> >
> > >
> > >
> > > Is there anything that I should consider more?
> > >
> > >
> > > 2016-08-19 4:47 GMT+09:00 Matthew Knepley <knepley at gmail.com>:
> > > On Thu, Aug 18, 2016 at 2:44 PM, 최경준 <kyungjun.choi92 at gmail.com>
> wrote:
> > > Is there a part that you considered this as finite-difference
> approximation?
> > > I thought I used matrix-free method with MatCreateSNESMF() function
> > >
> > > You did not tell the SNES to use a MF Jacobian, you just made a Mat
> object. This is why
> > > we encourage people to use the command line. Everything is setup
> correctly and in order.
> > > Why would you choose not to. This creates long rounds of email.
> > >
> > > Matt
> > >
> > > Also I used
> > > - call PCSetType(pc, PCNONE, ier) --> so the pc type shows 'none' at
> the log
> > >
> > >
> > > I didn't use any of command line options.
> > >
> > >
> > > Kyungjun
> > >
> > > 2016-08-19 4:27 GMT+09:00 Barry Smith <bsmith at mcs.anl.gov>:
> > >
> > > You can't use that Jacobian function SNESComputeJacobianDefault
> with matrix free, it tries to compute the matrix entries and stick them
> into the matrix. You can use MatMFFDComputeJacobian
> > >
> > > > On Aug 18, 2016, at 2:03 PM, 최경준 <kyungjun.choi92 at gmail.com> wrote:
> > > >
> > > > I got stuck at FormJacobian stage.
> > > >
> > > > - call SNESComputeJacobianDefault(snes, v, J, pJ, FormResidual,
> ier) --> J & pJ are same with A matrix-free matrix (input argument)
> > > >
> > > >
> > > >
> > > > with these kind of messages..
> > > >
> > > > [0]PETSC ERROR: No support for this operation for this object type
> > > > [0]PETSC ERROR: Mat type mffd
> > > >
> > > >
> > > >
> > > > Guess it's because I used A matrix-free matrix (which is mffd type)
> into pJ position.
> > > >
> > > > Is there any solution for this kind of situation?
> > > >
> > > >
> > > > 2016-08-19 2:05 GMT+09:00 Matthew Knepley <knepley at gmail.com>:
> > > > On Thu, Aug 18, 2016 at 12:04 PM, 최경준 <kyungjun.choi92 at gmail.com>
> wrote:
> > > > Then in order not to use preconditioner,
> > > >
> > > > is it ok if I just put A matrix-free matrix (made from
> MatCreateSNESMF()) into the place where preA should be?
> > > >
> > > > Yes, but again the solve will likely perform very poorly.
> > > >
> > > > Thanks,
> > > >
> > > > Matt
> > > >
> > > > The flow goes like this
> > > > - call SNESCreate
> > > > - call SNESSetFunction(snes, r, FormResidual, userctx, ier)
> > > > - call MatCreateSNESMF(snes, A, ier)
> > > > - call SNESSetJacobian(snes, A, A, FormJacobian, userctx, ier)
> > > > - call SNESSetFromOptions()
> > > >
> > > > - call SNESGetKSP(snes, ksp, ier)
> > > > - call KSPSetType(ksp, KSPGMRES, ier)
> > > > - call KSPGetPC(ksp, pc, ier)
> > > > - call PCSetType(pc, PCNONE, ier)
> > > > - call KSPGMRESSetRestart(ksp, 30, ier)
> > > >
> > > > - call SNESSolve()
> > > > .
> > > > .
> > > >
> > > >
> > > > and inside the FormJacobian routine
> > > > - call SNESComputeJacobian(snes, v, J, pJ, userctx, ier) --> J and
> pJ must be pointed with A and A.
> > > >
> > > >
> > > >
> > > > Thank you again,
> > > >
> > > > Kyungjun.
> > > >
> > > > 2016-08-19 1:44 GMT+09:00 Matthew Knepley <knepley at gmail.com>:
> > > > On Thu, Aug 18, 2016 at 11:42 AM, 최경준 <kyungjun.choi92 at gmail.com>
> wrote:
> > > > Thanks for your helpful answers.
> > > >
> > > > Here's another question...
> > > >
> > > > As I read some example PETSc codes, I noticed that there should be a
> preconditioning matrix (e.g. approx. jacobian matrix) when using
> MatCreateSNESMF().
> > > >
> > > > I mean,
> > > > after calling MatCreateSNESMF(snes, A, ier),
> > > > there should be another matrix preA(preconditioning matrix) to use
> SNESSetJacobian(snes, A, preA, FormJacobian, ctx, ier).
> > > >
> > > >
> > > > 1) Is there any way that I can use matrix-free method without making
> preconditioning matrix?
> > > >
> > > > Don't use a preconditioner. As you might expect, this does not often
> work out well.
> > > >
> > > > 2) I have a reference code, and the code adopts
> > > >
> > > > MatFDColoringCreate()
> > > > and finally uses
> > > > SNESComputeJacobianDefaultColor() at FormJacobian stage.
> > > >
> > > > But I can't see the inside of the fdcolor and I'm curious of this
> mechanism. Can you explain this very briefly or tell me an example code
> that I can refer to. ( I think none of PETSc example code is using
> fdcolor..)
> > > >
> > > > This is the default, so there is no need for all that code. We use
> naive graph 2-coloring. I think there might be a review article by Alex
> Pothen about that.
> > > >
> > > > Thanks,
> > > >
> > > > Matt
> > > >
> > > >
> > > > Best,
> > > >
> > > > Kyungjun.
> > > >
> > > > 2016-08-19 0:54 GMT+09:00 Matthew Knepley <knepley at gmail.com>:
> > > > On Thu, Aug 18, 2016 at 10:39 AM, 최경준 <kyungjun.choi92 at gmail.com>
> wrote:
> > > > 1) I wanna know the difference between applying option with command
> line and within source code.
> > > > From my experience, command line option helps set other default
> settings that I didn't applied, I guess.
> > > >
> > > > The command line arguments are applied to an object when
> *SetFromOptions() is called, so in this case
> > > > you want SNESSetFromOptions() on the solver. There should be no
> difference from using the API.
> > > >
> > > > 2) I made a matrix-free matrix with MatCreateSNESMF function, and
> every time I check my snes context with SNESView,
> > > >
> > > > Mat Object: 1 MPI processes
> > > > type: mffd
> > > > rows=11616, cols=11616
> > > > Matrix-free approximation:
> > > > err=1.49012e-08 (relative error in function evaluation)
> > > > The compute h routine has not yet been set
> > > >
> > > > at the end of line shows there's no routine for computing h value.
> > > > I used MatMFFDWPSetComputeNormU function, but it didn't work I think.
> > > > Is it ok if I leave the h value that way? Or should I have to set h
> computing routine?
> > > >
> > > > I am guessing you are calling the function on a different object
> from the one that is viewed here.
> > > > However, there will always be a default function for computing h.
> > > >
> > > > Thanks,
> > > >
> > > > Matt
> > > >
> > > > Kyungjun.
> > > >
> > > > 2016-08-18 23:18 GMT+09:00 Matthew Knepley <knepley at gmail.com>:
> > > > On Thu, Aug 18, 2016 at 8:35 AM, 최경준 <kyungjun.choi92 at gmail.com>
> wrote:
> > > > Hi, I'm trying to set my SNES matrix-free with Walker & Pernice way
> of computing h value.
> > > >
> > > > I found above command (MatSNESMFWPSetComputeNormU) but my fortran
> compiler couldn't fine any reference of that command.
> > > >
> > > > I checked Petsc changes log, but there weren't any mentions about
> that command.
> > > >
> > > > Should I have to include another specific header file?
> > > >
> > > > We have this function
> > > >
> > > > http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/
> MatMFFDWPSetComputeNormU.html
> > > >
> > > > but I would recommend using the command line option
> > > >
> > > > -mat_mffd_compute_normu
> > > >
> > > > Thanks,
> > > >
> > > > Matt
> > > >
> > > > Thank you always.
> > > >
> > > >
> > > >
> > > > --
> > > > 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
> > > >
> > > >
> > > >
> > > >
> > > > --
> > > > 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
> > > >
> > > >
> > > >
> > > >
> > > > --
> > > > 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
> > > >
> > > >
> > > >
> > > >
> > > > --
> > > > 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
> > > >
> > >
> > >
> > >
> > >
> > >
> > > --
> > > 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
> > >
> >
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20161006/f716392c/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 22943 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20161006/f716392c/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 222628 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20161006/f716392c/attachment-0003.png>
More information about the petsc-users
mailing list