[petsc-users] Question about using MatSNESMFWPSetComputeNormU
최경준
kyungjun.choi92 at gmail.com
Fri Aug 19 00:04:44 CDT 2016
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/20160819/d3eac1dc/attachment-0001.html>
More information about the petsc-users
mailing list