[petsc-users] Question about using MatSNESMFWPSetComputeNormU
Matthew Knepley
knepley at gmail.com
Fri Aug 19 06:47:56 CDT 2016
On Fri, 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)*
>
Remove these two lines
> *call MatCreateSNESMF(Mixt%snes, Mixt%A, ier)*
>
> *call SNESSetJacobian(Mixt%snes, Mixt%A, Mixt%A, MatMFFDComputeJacobian,
> Collect, ier)*
>
> *call SNESSetFromOptions(Mixt%snes, ier)*
>
And also all this stuff
> *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)*
>
until here. Then give
-snes_mf -pc_type none -snes_view -snes_monitor
-ksp_monitor_true_residual -ksp_converged_reason -snes_convergred_reason
and you have what you want.
Matt
> *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
>> >
>>
>>
>
--
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/ac32e9ef/attachment-0001.html>
More information about the petsc-users
mailing list