> On Jul 14, 2021, at 5:01 PM, Tang, Qi <[email protected]> wrote: > > Thanks a lot for the explanation, Matt and Stefano. That helps a lot. > > Just to confirm, the comment in src/ts/impls/implicit/theta/theta.c seems to > indicates TS solves U_{n+1} in its SNES/KSP solve, but it actually solves > the update dU_n in U_{n+1} = U_n - lambda*dU_n in the solve. Right?
The SNES object solves the nonlinear equations as written in the comment of TSTHETA. F[t0+Theta*dt, U, (U-U0)*shift] = 0 In case SNES is of type SNESLS (Newton), then the linearized equations are solved. The linear system matrix is the one provided by the IJacobian function J = dF/dU + shift dF/dUdot If it is SNESKSPONLY ( as it should be for TS_LINEAR), then only one step is taken and lambda = 1. > > It actually makes a lot sense, because KSPSolve in TSSolve reports it uses > zero initial guess. So if what I said is true, that effectively means it uses > U0 as the initial guess. > > Qi > >> On Jul 14, 2021, at 2:56 AM, Matthew Knepley <[email protected] >> <mailto:[email protected]>> wrote: >> >> On Wed, Jul 14, 2021 at 4:43 AM Stefano Zampini <[email protected] >> <mailto:[email protected]>> wrote: >> Qi >> >> Backward Euler is a special case of Theta methods in PETSc (Theta=1). In >> src/ts/impls/implicit/theta/theta.c on top of SNESTSFormFunction_Theta you >> have some explanation of what is solved for at each time step (see below). >> SNES then solves for the Newton update dy_n and the next Newton iterate is >> computed as x_{n+1} = x_{n} - lambda * dy_n. Hope this helps. >> >> In other words, you should be able to match the initial residual to >> >> F(t + dt, 0, -Un / dt) >> >> for your IFunction. However, it is really not normal to use U = 0. The >> default is to use U = U0 >> as the initial guess I think. >> >> Thanks, >> >> Matt >> >> /* >> This defines the nonlinear equation that is to be solved with SNES >> G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 >> >> Note that U here is the stage argument. This means that U = U_{n+1} only >> if endpoint = true, >> otherwise U = theta U_{n+1} + (1 - theta) U0, which for the case of >> implicit midpoint is >> U = (U_{n+1} + U0)/2 >> */ >> static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts) >> >> >>> On Jul 14, 2021, at 6:12 AM, Tang, Qi <[email protected] >>> <mailto:[email protected]>> wrote: >>> >>> Hi, >>> >>> During the process to experiment the suggestion Matt made, we ran into some >>> questions regarding to TSSolve vs KSPSolve. We got different initial >>> unpreconditioned residual using two solvers. Let’s say we solve the problem >>> with backward Euler and there is no rhs. We guess TSSolve solves >>> (U^{n+1}-U^n)/dt = A U^{n+1}. >>> (We only provides IJacobian in this case and turn on TS_LINEAR.) >>> So we guess the initial unpreconditioned residual would be ||U^n/dt||_2, >>> which seems different from the residual we got from a backward Euler >>> stepping we implemented by ourself through KSPSolve. >>> >>> Do we have some misunderstanding on TSSolve? >>> >>> Thanks, >>> Qi >>> T5@LANL >>> >>> >>> >>>> On Jul 7, 2021, at 3:54 PM, Matthew Knepley <[email protected] >>>> <mailto:[email protected]>> wrote: >>>> >>>> On Wed, Jul 7, 2021 at 2:33 PM Jorti, Zakariae <[email protected] >>>> <mailto:[email protected]>> wrote: >>>> Hi Matt, >>>> >>>> >>>> >>>> Thanks for your quick reply. >>>> >>>> I have not completely understood your suggestion, could you please >>>> elaborate a bit more? >>>> >>>> For your convenience, here is how I am proceeding for the moment in my >>>> code: >>>> >>>> >>>> >>>> TSGetKSP(ts,&ksp); >>>> >>>> KSPGetPC(ksp,&pc); >>>> >>>> PCSetType(pc,PCFIELDSPLIT); >>>> >>>> PCFieldSplitSetDetectSaddlePoint(pc,PETSC_TRUE); >>>> >>>> PCSetUp(pc); >>>> >>>> PCFieldSplitGetSubKSP(pc, &n, &subksp); >>>> >>>> KSPGetPC(subksp[1], &(subpc[1])); >>>> >>>> I do not like the two lines above. We should not have to do this. >>>> KSPSetOperators(subksp[1],T,T); >>>> >>>> In the above line, I want you to use a separate preconditioning matrix M, >>>> instead of T. That way, it will provide >>>> the preconditioning matrix for your Schur complement problem. >>>> >>>> Thanks, >>>> >>>> Matt >>>> KSPSetUp(subksp[1]); >>>> >>>> PetscFree(subksp); >>>> >>>> TSSolve(ts,X); >>>> >>>> >>>> >>>> Thank you. >>>> >>>> Best, >>>> >>>> >>>> >>>> Zakariae >>>> >>>> From: Matthew Knepley <[email protected] <mailto:[email protected]>> >>>> Sent: Wednesday, July 7, 2021 12:11:10 PM >>>> To: Jorti, Zakariae >>>> Cc: [email protected] <mailto:[email protected]>; Tang, Qi; >>>> Tang, Xianzhu >>>> Subject: [EXTERNAL] Re: [petsc-users] Problem with PCFIELDSPLIT >>>> >>>> On Wed, Jul 7, 2021 at 1:51 PM Jorti, Zakariae via petsc-users >>>> <[email protected] <mailto:[email protected]>> wrote: >>>> Hi, >>>> >>>> >>>> >>>> I am trying to build a PCFIELDSPLIT preconditioner for a matrix >>>> >>>> J = [A00 A01] >>>> >>>> [A10 A11] >>>> >>>> that has the following shape: >>>> >>>> >>>> >>>> M_{user}^{-1} = [I -ksp(A00) A01] [ksp(A00) 0] [I >>>> 0] >>>> >>>> [0 I] [0 >>>> ksp(T)] [-A10 ksp(A00) I ] >>>> >>>> >>>> >>>> where T is a user-defined Schur complement approximation that replaces the >>>> true Schur complement S:= A11 - A10 ksp(A00) A01. >>>> >>>> >>>> >>>> I am trying to do something similar to this example (lines 41--45 and >>>> 116--121): >>>> https://www.mcs.anl.gov/petsc/petsc-current/src/snes/tutorials/ex70.c.html >>>> <https://urldefense.com/v3/__https://www.mcs.anl.gov/petsc/petsc-current/src/snes/tutorials/ex70.c.html__;!!HXCxUKc!hoEfgnaraTfQoSgAiplsc6GJ_HuPXN88m5AJVy1gb7WVMNkGENDnJ3zToOGlhw$> >>>> >>>> The problem I have is that I manage to replace S with T on a separate >>>> single linear system but not for the linear systems generated by my >>>> time-dependent PDE. Even if I set the preconditioner M_{user}^{-1} >>>> correctly, the T matrix gets replaced by S in the preconditioner once I >>>> call TSSolve. >>>> >>>> Do you have any suggestions how to fix this knowing that the matrix J does >>>> not change over time? >>>> >>>> >>>> I don't like how it is done in that example for this very reason. >>>> >>>> When I want to use a custom preconditioning matrix for the Schur >>>> complement, I always give a preconditioning matrix M to the outer solve. >>>> Then PCFIELDSPLIT automatically pulls the correct block from M, (1,1) for >>>> the Schur complement, for that preconditioning matrix without >>>> extra code. Can you do this? >>>> >>>> Thanks, >>>> >>>> Matt >>>> Many thanks. >>>> >>>> >>>> >>>> Best regards, >>>> >>>> >>>> >>>> Zakariae >>>> >>>> >>>> >>>> >>>> -- >>>> 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/ >>>> <https://urldefense.com/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!HXCxUKc!hoEfgnaraTfQoSgAiplsc6GJ_HuPXN88m5AJVy1gb7WVMNkGENDnJ3wB3dcMFw$> >>>> >>>> >>>> -- >>>> 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/ >>>> <https://urldefense.com/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!HXCxUKc!hoEfgnaraTfQoSgAiplsc6GJ_HuPXN88m5AJVy1gb7WVMNkGENDnJ3wB3dcMFw$> >> >> >> >> -- >> 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/ >> <https://urldefense.com/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!HXCxUKc!msQrz7__TrpOmaTVhvY1yLAlDQXNJ5jcYVAxF4lcpyLrZqt2lFe22bkbuJMizA$>
