<div dir="ltr">Dear Barry and Matt.<div><br></div><div>Thank you very much for helping me up all night. (in my time)</div><div><br></div><div>And sorry for not asking with sufficient source code condition or my circumstances. (also with poor English.)</div><div><br></div><div><br></div><div>I just want to make sure that the options of my code is well applied.</div><div><br></div><div>I'm trying to use GMRES with matrix-free method. I'd like to solve 2-D euler equation without preconditioning matrix, for now.</div><div><br></div><div><br></div><div>1) I'm still curious whether my snes context is using MF jacobian. ( just like -snes_mf command line option)</div><div><br></div><div>2) And mind if I ask you that whether I applied petsc functions properly?</div><div><br></div><div>I'll check out ex5 for applying command line options.</div><div><br></div><div><br></div><div>I'll attach my petsc flow code and option log by SNESView() below.</div><div>--------------------------------------------------------------------------------------------------------------------<br></div><div><b>- petsc flow code</b></div><div>--------------------------------------------------------------------------------------------------------------------<br></div><div><br></div><div>ndof = Mixt%nCVar * Grid%nCell<br></div><div><br></div><div><b>call VecCreateMPIWIthArray(PETSC_COMM_WORLD, Mixt%nCVar, ndof, PETSC_DECIDE, Mixt%cv, Mixt%x, ier)</b></div><div><b>call VecDuplicate(Mixt%x, Mixt%r, ier)</b></div><div><b>call VecSet(Mixt%r, zero, ier)</b></div><div><b><br></b></div><div><b>call SNESCreate(PETSC_COMM_WORLD, Mixt%snes, ier)</b></div><div><b>call SNESSetFunction(Mixt%snes, Mixt%r, FormPetscResidual, Collect, ier)</b></div><div><b>call MatCreateSNESMF(Mixt%snes, Mixt%A, ier)</b></div><div><b><br></b></div><div><b>call SNESSetJacobian(Mixt%snes, Mixt%A, Mixt%A, MatMFFDComputeJacobian, Collect, ier)</b></div><div><b>call SNESSetFromOptions(Mixt%snes, ier)</b></div><div><b><br></b></div><div><b>call SNESGetKSP(Mixt%snes, ksp, ier)</b></div><div><b>call KSPSetType(ksp, KSPGMRES, ier)</b></div><div><b>call KSPGetPC(ksp, pc, ier)</b></div><div><b>call PCSetType(pc, PCNONE, ier)</b></div><div><b>call KSPSetInitialGuessNonzero(ksp, PETSC_TRUE, ier)</b></div><div><b>call KSPGMRESSetRestart(ksp, 30, ier)</b></div><div><b>call KSPGMRESSetPreAllocation(ksp, ier)</b></div><div><b><br></b></div><div><b><br></b></div><div><div><div><b>call SNESSetFunction(Mixt%snes, Mixt%r, FormPetscResidual, Collect, ier)</b></div><div><b>call SNESSetJacobian(Mixt%snes, Mixt%A, Mixt%A, MatMFFDComputeJacobian, Collect, ier)</b></div></div><div><b><br></b></div><div><b>call SNESSolve(Mixt%snes, PETSC_NULL_OBJECT, Mixt%x, ier)</b></div></div><div><br></div><div><b>stop</b> ( for temporary )</div><div><br></div><div><br></div><div>--------------------------------------------------------------------------------------------------------------------</div><div><b>subroutine FormPetscResidual(snes, x, f, Collect, ier)</b></div><div><div> type(t_Collect), intent(inout) :: Collect</div><div><br></div><div> SNES :: snes</div><div> Vec :: x, f</div><div> integer :: ier, counter, iCell, iVar, temp</div><div> integer :: ndof</div><div> real(8), allocatable :: CVar(:,:)</div><div> real(8), allocatable :: PVar(:,:)</div><div> PetscScalar, pointer :: xx_v(:)</div><div> PetscScalar, pointer :: ff_v(:) </div><div> </div><div> ! Set degree of freedom of this system.</div><div> ndof = Collect%pMixt%nCVar * Collect%pGrid%nCell</div><div><br></div><div> ! Backup the original values for cv to local array CVar</div><div> allocate( CVar(0:Collect%pMixt%nCVar-1, Collect%pGrid%nCell) )</div><div> allocate( PVar(0:Collect%pMixt%nPVar-1, Collect%pGrid%nCell) )</div><div> allocate( xx_v(1:ndof) )</div><div> allocate( ff_v(1:ndof) )</div><div> xx_v(:) = 0d0</div><div> ff_v(:) = 0d0</div><div><br></div><div> ! Backup the original values for cv and pv</div><div> do iCell = 1, Collect%pGrid%nCell</div><div> do iVar = 0, Collect%pMixt%nCVar-1</div><div> CVar(iVar,iCell) = Collect%pMixt%cv(iVar,iCell)</div><div> PVar(iVar,iCell) = Collect%pMixt%pv(iVar,iCell)</div><div> end do</div><div> end do</div><div><br></div><div> ! Copy the input argument vector x to array value xx_v</div><div> call VecGetArrayReadF90(x, xx_v, ier)</div><div> call VecGetArrayF90(f, ff_v, ier)</div><div><br></div><div> ! Compute copy the given vector into Mixt%cv and check for validity</div><div> counter = 0</div><div> do iCell = 1, Collect%pGrid%nCell</div><div> do iVar = 0, Collect%pMixt%nCVar-1</div><div> counter = counter + 1</div><div> Collect%pMixt%cv(iVar,iCell) = xx_v(counter)</div><div> end do</div><div> end do</div><div> </div><div> ! Update primitive variables with input x vector to compute residual</div><div> call PostProcessing(Collect%pMixt,Collect%pGrid,Collect%pConf)</div><div><br></div><div><br></div><div> ! Compute the residual</div><div> call ComputeResidual(Collect%pMixt,Collect%pGrid,Collect%pConf) --> where update residual of cell</div><div> </div><div> ! Copy the residual array into the PETSc vector</div><div> counter = 0</div><div> do iCell = 1, Collect%pGrid%nCell</div><div> do iVar = 0, Collect%pMixt%nCVar-1</div><div> counter = counter + 1</div><div><br></div><div> ff_v(counter) = Collect%pMixt%Residual(iVar,iCell) + Collect%pGrid%vol(iCell)/Collect%pMixt%TimeStep(iCell)*( Collect%pMixt%cv(iVar,iCell) - CVar(iVar,iCell) )</div><div> end do</div><div> end do</div><div><br></div><div> ! Restore conservative variables</div><div> do iCell = 1, Collect%pGrid%nCell</div><div> do iVar = 0, Collect%pMixt%nCVar-1</div><div> Collect%pMixt%cv(iVar,iCell) = CVar(iVar,iCell)</div><div> Collect%pMixt%pv(iVar,iCell) = PVar(iVar,iCell)</div><div> end do</div><div> end do</div><div><br></div><div> call VecRestoreArrayReadF90(x, xx_v, ier)</div><div> call VecRestoreArrayF90(f, ff_v, ier)</div><div><br></div><div> deallocate(CVar)</div><div> deallocate(PVar)</div></div><div>--------------------------------------------------------------------------------------------------------------------<br></div><div><br></div><div><br></div><div>--------------------------------------------------------------------------------------------------------------------<br></div><div><b>- option log</b></div><div>--------------------------------------------------------------------------------------------------------------------</div><div><div>SNES Object: 1 MPI processes</div><div> type: newtonls</div><div> SNES has not been set up so information may be incomplete</div><div> maximum iterations=1, maximum function evaluations=10000</div><div> tolerances: relative=1e-08, absolute=1e-32, solution=1e-08</div><div> total number of linear solver iterations=0</div><div> total number of function evaluations=0</div><div> norm schedule ALWAYS</div><div> SNESLineSearch Object: 1 MPI processes</div><div> type: bt</div><div> interpolation: cubic</div><div> alpha=1.000000e-04</div><div> maxstep=1.000000e+08, minlambda=1.000000e-12</div><div> tolerances: relative=1.000000e-08, absolute=1.000000e-15, lambda=1.000000e-08</div><div> maximum iterations=40</div><div> KSP Object: 1 MPI processes</div><div> type: gmres</div><div> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement</div><div> GMRES: happy breakdown tolerance 1e-30</div><div> maximum iterations=10000</div><div> tolerances: relative=0.001, absolute=1e-50, divergence=10000.</div><div> left preconditioning</div><div> using nonzero initial guess</div><div> using DEFAULT norm type for convergence test</div><div> PC Object: 1 MPI processes</div><div> type: none</div><div> PC has not been set up so information may be incomplete</div><div> linear system matrix = precond matrix:</div><div> Mat Object: 1 MPI processes</div><div> type: mffd</div><div> rows=11616, cols=11616</div><div> Matrix-free approximation:</div><div> err=1.49012e-08 (relative error in function evaluation)</div><div> The compute h routine has not yet been set</div></div><div><br></div><div><br></div><div>Sincerely,</div><div><br></div><div>Kyungjun</div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">2016-08-19 13:00 GMT+09:00 Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class=""><br>
> On Aug 18, 2016, at 10:28 PM, 최경준 <<a href="mailto:kyungjun.choi92@gmail.com">kyungjun.choi92@gmail.com</a>> wrote:<br>
><br>
> Dear Matt.<br>
><br>
> I didn't use the command line options because it looked not working.<br>
><br>
> I called SNESSetFromOptions(snes, ier) in my source code,<br>
><br>
> but options like -snes_mf or -snes_monitor doesn't look working.<br>
<br>
</span> "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".<br>
<span class="HOEnZb"><font color="#888888"><br>
Barry<br>
</font></span><div class="HOEnZb"><div class="h5"><br>
<br>
><br>
><br>
> Is there anything that I should consider more?<br>
><br>
><br>
> 2016-08-19 4:47 GMT+09:00 Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>>:<br>
> On Thu, Aug 18, 2016 at 2:44 PM, 최경준 <<a href="mailto:kyungjun.choi92@gmail.com">kyungjun.choi92@gmail.com</a>> wrote:<br>
> Is there a part that you considered this as finite-difference approximation?<br>
> I thought I used matrix-free method with MatCreateSNESMF() function<br>
><br>
> You did not tell the SNES to use a MF Jacobian, you just made a Mat object. This is why<br>
> we encourage people to use the command line. Everything is setup correctly and in order.<br>
> Why would you choose not to. This creates long rounds of email.<br>
><br>
> Matt<br>
><br>
> Also I used<br>
> - call PCSetType(pc, PCNONE, ier) --> so the pc type shows 'none' at the log<br>
><br>
><br>
> I didn't use any of command line options.<br>
><br>
><br>
> Kyungjun<br>
><br>
> 2016-08-19 4:27 GMT+09:00 Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>>:<br>
><br>
> 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<br>
><br>
> > On Aug 18, 2016, at 2:03 PM, 최경준 <<a href="mailto:kyungjun.choi92@gmail.com">kyungjun.choi92@gmail.com</a>> wrote:<br>
> ><br>
> > I got stuck at FormJacobian stage.<br>
> ><br>
> > - call SNESComputeJacobianDefault(<wbr>snes, v, J, pJ, FormResidual, ier) --> J & pJ are same with A matrix-free matrix (input argument)<br>
> ><br>
> ><br>
> ><br>
> > with these kind of messages..<br>
> ><br>
> > [0]PETSC ERROR: No support for this operation for this object type<br>
> > [0]PETSC ERROR: Mat type mffd<br>
> ><br>
> ><br>
> ><br>
> > Guess it's because I used A matrix-free matrix (which is mffd type) into pJ position.<br>
> ><br>
> > Is there any solution for this kind of situation?<br>
> ><br>
> ><br>
> > 2016-08-19 2:05 GMT+09:00 Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>>:<br>
> > On Thu, Aug 18, 2016 at 12:04 PM, 최경준 <<a href="mailto:kyungjun.choi92@gmail.com">kyungjun.choi92@gmail.com</a>> wrote:<br>
> > Then in order not to use preconditioner,<br>
> ><br>
> > is it ok if I just put A matrix-free matrix (made from MatCreateSNESMF()) into the place where preA should be?<br>
> ><br>
> > Yes, but again the solve will likely perform very poorly.<br>
> ><br>
> > Thanks,<br>
> ><br>
> > Matt<br>
> ><br>
> > The flow goes like this<br>
> > - call SNESCreate<br>
> > - call SNESSetFunction(snes, r, FormResidual, userctx, ier)<br>
> > - call MatCreateSNESMF(snes, A, ier)<br>
> > - call SNESSetJacobian(snes, A, A, FormJacobian, userctx, ier)<br>
> > - call SNESSetFromOptions()<br>
> ><br>
> > - call SNESGetKSP(snes, ksp, ier)<br>
> > - call KSPSetType(ksp, KSPGMRES, ier)<br>
> > - call KSPGetPC(ksp, pc, ier)<br>
> > - call PCSetType(pc, PCNONE, ier)<br>
> > - call KSPGMRESSetRestart(ksp, 30, ier)<br>
> ><br>
> > - call SNESSolve()<br>
> > .<br>
> > .<br>
> ><br>
> ><br>
> > and inside the FormJacobian routine<br>
> > - call SNESComputeJacobian(snes, v, J, pJ, userctx, ier) --> J and pJ must be pointed with A and A.<br>
> ><br>
> ><br>
> ><br>
> > Thank you again,<br>
> ><br>
> > Kyungjun.<br>
> ><br>
> > 2016-08-19 1:44 GMT+09:00 Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>>:<br>
> > On Thu, Aug 18, 2016 at 11:42 AM, 최경준 <<a href="mailto:kyungjun.choi92@gmail.com">kyungjun.choi92@gmail.com</a>> wrote:<br>
> > Thanks for your helpful answers.<br>
> ><br>
> > Here's another question...<br>
> ><br>
> > As I read some example PETSc codes, I noticed that there should be a preconditioning matrix (e.g. approx. jacobian matrix) when using MatCreateSNESMF().<br>
> ><br>
> > I mean,<br>
> > after calling MatCreateSNESMF(snes, A, ier),<br>
> > there should be another matrix preA(preconditioning matrix) to use SNESSetJacobian(snes, A, preA, FormJacobian, ctx, ier).<br>
> ><br>
> ><br>
> > 1) Is there any way that I can use matrix-free method without making preconditioning matrix?<br>
> ><br>
> > Don't use a preconditioner. As you might expect, this does not often work out well.<br>
> ><br>
> > 2) I have a reference code, and the code adopts<br>
> ><br>
> > MatFDColoringCreate()<br>
> > and finally uses<br>
> > SNESComputeJacobianDefaultColo<wbr>r() at FormJacobian stage.<br>
> ><br>
> > 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..)<br>
> ><br>
> > 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.<br>
> ><br>
> > Thanks,<br>
> ><br>
> > Matt<br>
> ><br>
> ><br>
> > Best,<br>
> ><br>
> > Kyungjun.<br>
> ><br>
> > 2016-08-19 0:54 GMT+09:00 Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>>:<br>
> > On Thu, Aug 18, 2016 at 10:39 AM, 최경준 <<a href="mailto:kyungjun.choi92@gmail.com">kyungjun.choi92@gmail.com</a>> wrote:<br>
> > 1) I wanna know the difference between applying option with command line and within source code.<br>
> > From my experience, command line option helps set other default settings that I didn't applied, I guess.<br>
> ><br>
> > The command line arguments are applied to an object when *SetFromOptions() is called, so in this case<br>
> > you want SNESSetFromOptions() on the solver. There should be no difference from using the API.<br>
> ><br>
> > 2) I made a matrix-free matrix with MatCreateSNESMF function, and every time I check my snes context with SNESView,<br>
> ><br>
> > Mat Object: 1 MPI processes<br>
> > type: mffd<br>
> > rows=11616, cols=11616<br>
> > Matrix-free approximation:<br>
> > err=1.49012e-08 (relative error in function evaluation)<br>
> > The compute h routine has not yet been set<br>
> ><br>
> > at the end of line shows there's no routine for computing h value.<br>
> > I used MatMFFDWPSetComputeNormU function, but it didn't work I think.<br>
> > Is it ok if I leave the h value that way? Or should I have to set h computing routine?<br>
> ><br>
> > I am guessing you are calling the function on a different object from the one that is viewed here.<br>
> > However, there will always be a default function for computing h.<br>
> ><br>
> > Thanks,<br>
> ><br>
> > Matt<br>
> ><br>
> > Kyungjun.<br>
> ><br>
> > 2016-08-18 23:18 GMT+09:00 Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>>:<br>
> > On Thu, Aug 18, 2016 at 8:35 AM, 최경준 <<a href="mailto:kyungjun.choi92@gmail.com">kyungjun.choi92@gmail.com</a>> wrote:<br>
> > Hi, I'm trying to set my SNES matrix-free with Walker & Pernice way of computing h value.<br>
> ><br>
> > I found above command (MatSNESMFWPSetComputeNormU) but my fortran compiler couldn't fine any reference of that command.<br>
> ><br>
> > I checked Petsc changes log, but there weren't any mentions about that command.<br>
> ><br>
> > Should I have to include another specific header file?<br>
> ><br>
> > We have this function<br>
> ><br>
> > <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMFFDWPSetComputeNormU.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-current/docs/<wbr>manualpages/Mat/<wbr>MatMFFDWPSetComputeNormU.html</a><br>
> ><br>
> > but I would recommend using the command line option<br>
> ><br>
> > -mat_mffd_compute_normu<br>
> ><br>
> > Thanks,<br>
> ><br>
> > Matt<br>
> ><br>
> > Thank you always.<br>
> ><br>
> ><br>
> ><br>
> > --<br>
> > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
> > -- Norbert Wiener<br>
> ><br>
> ><br>
> ><br>
> ><br>
> > --<br>
> > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
> > -- Norbert Wiener<br>
> ><br>
> ><br>
> ><br>
> ><br>
> > --<br>
> > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
> > -- Norbert Wiener<br>
> ><br>
> ><br>
> ><br>
> ><br>
> > --<br>
> > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
> > -- Norbert Wiener<br>
> ><br>
><br>
><br>
><br>
><br>
><br>
> --<br>
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
> -- Norbert Wiener<br>
><br>
<br>
</div></div></blockquote></div><br></div>