[petsc-users] Implementing TS routine
Manuel Valera
mvalera-w at sdsu.edu
Tue Jul 9 18:30:34 CDT 2019
Ok thanks everyone for your replies, Matt thanks for suggesting TStep, it
solved my problem for reusing the TS objects.
The TS routine are working as intended, and i was able to implement the TS
for the velocities as well, i dropped the use of the natural vector and i
am using just the DM global vectors as function and input for the TS RHS
function.
I'll write back if i have further questions,
Thanks so much,
On Tue, Jul 9, 2019 at 1:32 PM Matthew Knepley <knepley at gmail.com> wrote:
> On Tue, Jul 9, 2019 at 2:39 PM Manuel Valera via petsc-users <
> petsc-users at mcs.anl.gov> wrote:
>
>> On Tue, Jul 9, 2019 at 11:27 AM Smith, Barry F. <bsmith at mcs.anl.gov>
>> wrote:
>>
>>> > On Jul 8, 2019, at 6:53 PM, Manuel Valera via petsc-users <
>>> petsc-users at mcs.anl.gov> wrote:
>>> >
>>> > Hi Zhang,
>>> >
>>> > Thanks to your help i have implemented the TS routine for my
>>> temperature DMDA array, and it works correctly including in MPI.
>>> >
>>> > The breakthrough for me was realizing that TS steps and max time steps
>>> are all set up once, and used/advanced all at once when invoked with
>>> TSSolution. So what i did was to add all of the TSCreate / initialization
>>> into the main loop so i am basically creating and destroying the TS objects
>>> every time step.
>>>
>>> Why? You can create the TS objects once outside and then inside
>>> advance them one step at a time inside the loop if you want.
>>>
>>
>> For some reason this is not possible, if you could provide any ideas let
>> me know, the situation is this:
>>
>> This block is now inside the main loop:
>> call TSCreate(PETSC_COMM_WORLD,ts,ierr)
>> call TSSetProblemType(ts,TS_NONLINEAR,ierr)
>> call TSSetSolution(ts,LocTemperature,ierr)
>> call
>> TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,PETSC_NULL_VEC,ierr)
>> call TSSetType(ts,TSSSP,ierr)
>> call TSSetTimeStep(ts,dt,ierr)
>> call TSSetDM(ts,daScalars,ierr)
>> call TSSetMaxTime(ts,dt,ierr)
>> call
>> TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP,ierr)
>> call TSSetFromOptions(ts,ierr)
>> call TSSetUp(ts,ierr)
>> call TSSolve(ts,LocTemperature,ierr)
>> call TSDestroy(ts,ierr)
>>
>> Where i create, set everything up and solve the TS problem for a single
>> time step.
>>
>
> I admit that I have not completely followed the discussion. I will try and
> ask a few questions so that I can understand.
>
> First, you can use TSStep(), instead of TSSolve(), in order to run a
> single time step.
>
>
>> If i take TSCreate and TSDestroy out of this loop, the TS will advance
>> all the dt's at time zero of the program (first iteration) and then
>> continue executing the rest without taking this into account. The only way
>> i have found for TS to be executed every time in the way i need
>>
>
> Is it easy to say what you need to do between steps?
>
>
>> it is creating all objects and executing them at this moment. Any idea on
>> how i can control this from out of the loop would be appreciated.
>>
>> Keep in mind that i want to use only the TS routine in this step, as i
>> mentioned every example seems to include the dynamics of the problem inside
>> the TS, and then it makes sense that we could set a specific totaltime and
>> advance it all at TSSolution, what i am trying to do is quite different, as
>> all of the dynamics are to be preserved outside of the TS, i only really
>> need to make use of the function that takes u and f(t,u) and gives out
>> u(t+1) with a specific time stepping scheme.
>>
>>>
>>> > I hope the overhead time is not to bad when i scale to bigger problems.
>>>
>>> I would say it is only extremely strange conditions would you need to
>>> recreate the TS for each timestep. Better to now figure out how to reuse
>>> them.
>>>
>>
>> Agreed and i would like to know how to reuse them but the behavior seems
>> to be against what i want to do.
>>
>>
>> >
>>> > Next step would be to apply the same to advance the velocities, but
>>> problem now is that they are in a 4D dmda array of 3 degrees of freedom,
>>> any suggestions on how to implement this? does TS support multiple degrees
>>> of freedom arrays?
>>>
>>> TS doesn't care how you access the vector entries in your function
>>> evaluations. All TS every sees is a Vec.
>>>
>>> I am scared that you said 4D dmda array of 3 degrees of freedom.
>>> PETSc DMDA only goes up to 3 dimensions. So if you truely need 4 dim plus 3
>>> dof per point you need to "fake it" but putting the last dimension into the
>>> dof slot. In other words use a 3d DMDA and have for dof 3*n where n is the
>>> size of the 4th dimension. The code to access the values will now be a
>>> little cumbersome since you can only access the array as 4d and the 4th
>>> dimension contains both the 3 and the last dimension but it is not too bad.
>>>
>>> But if you mean 3 physical dimensions plus 3 dof then everything is
>>> fine just use the 3d DMDA.
>>>
>>
>> Sorry i didn't speak clearly, i really meant 3 physical dimensions plus 3
>> DOF, i am already using the 3D DMDAs here. My velocity array looks like
>> vel(0:3,i,j,k) so 3 spatial dimensions and 3 degrees of freedom one for
>> each velocity direction, so vel(0,...) is u, vel(1,..) is v and vel(2,..)
>> is w. The problem comes again when i try to use the TS to solve for this
>> kind of array. Do i need a 1 dimensional vector for the f(t,u) ?
>>
>
> There should be no problem with multiple DOF per point in your domain. You
> just take this into account
> in your residual (and Jacobian) function. There are many examples in TS
> with multiple DOFs. Maybe there
> is something I am missing here.
>
> Thanks,
>
> Matt
>
>
>> In the case of the temperature, i have created a natural vector from the
>> 3D RHS array, and use that as a f(t,u), do i need to do the same for the
>> velocities, or can i provide a 3D array with 3 DOF as either the f(t,u) and
>> the u? If the latter is the case, i could create a new 3D DMDA to hold the
>> RHS functions of each of the velocity, in the same structure the velocities
>> are. If the case if the former i would need to create new natural vectors
>> for each of the RHS, but also copy the velocities into a single DOF DMDA ?
>>
>> Thanks, i think i am a bit at a loss here,
>>
>> Regards,
>>
>>
>>
>>
>>
>>
>>
>>>
>>> Barry
>>>
>>>
>>>
>>> >
>>> > Thanks,
>>> >
>>> > On Thu, Jul 4, 2019 at 9:18 PM Zhang, Hong <hongzhang at anl.gov> wrote:
>>> >
>>> >
>>> >> On Jul 3, 2019, at 3:10 PM, Manuel Valera <mvalera-w at sdsu.edu> wrote:
>>> >>
>>> >> Thanks Zhang for your answer,
>>> >>
>>> >> I ended up getting a compiling and running TS routine... that does
>>> not give me the answers, and i am not sure what i am doing wrong,
>>> >>
>>> >> My TS code so far looks like this:
>>> >>
>>> >> (...initialization...)
>>> >>
>>> >> call TSCreate(PETSC_COMM_WORLD,ts,ierr)
>>> >> call TSSetProblemType(ts,TS_NONLINEAR,ierr)
>>> >> call TSSetSolution(ts,gTemperature,ierr)
>>> >> call
>>> TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,PETSC_NULL_VEC,ierr)
>>> >
>>> > The second last argument should be the user context. If you pass NULL,
>>> the ctx variable in your FormRHSFunction will be NULL.
>>> >
>>> >> call TSSetType(ts,TSSSP,ierr)
>>> >> call TSSetTimeStep(ts,dt,ierr)
>>> >> call TSSetDM(ts,daScalars,ierr)
>>> >> call TSSetMaxTime(ts,TotalTime,ierr)
>>> >> call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)
>>> >> call TSSetFromOptions(ts,ierr)
>>> >> call TSSetUp(ts,ierr)
>>> >>
>>> >> (... then inside main calculations loop i have...)
>>> >>
>>> >> call TSSolve(ts,gTemperature,ierr)
>>> >>
>>> >> (...and my RHSFunction looks like this...)
>>> >>
>>> >> subroutine FormRHSFunction(ts,t,t_loc,rhs_loc,ctx,ierr)
>>> >>
>>> >> use petscvec
>>> >> use petscts
>>> >> use petscdmda
>>> >> use petscdm
>>> >>
>>> >> USE utils
>>> >> USE dmdaobjs
>>> >> USE probsize
>>> >> USE modelparams
>>> >>
>>> >> implicit none
>>> >>
>>> >> TS :: ts
>>> >> PetscReal :: t
>>> >> Vec :: ctx,t_loc,rhs_loc
>>> >> PetscErrorCode :: ierr
>>> >>
>>> >>
>>> >>
>>> >> call TSGetDM(ts,daScalars,ierr)
>>> >>
>>> >> call DMGetLocalVector(daScalars,t_loc,ierr);
>>> >
>>> > t_loc is the input, it should be not modified.
>>> >
>>> >> call
>>> DMGlobalToLocalBegin(daScalars,gTemperature,INSERT_VALUES,t_loc,ierr);CHKERRQ(ierr)
>>> >> call
>>> DMGlobalToLocalEnd(daScalars,gTemperature,INSERT_VALUES,t_loc,ierr);CHKERRQ(ierr)
>>> >
>>> > FormRHSFunction is supposed to implement rhs_loc = f(t,t_loc). So you
>>> want to scatter ghost points for t_loc, not gTemperature.
>>> >
>>> >> call DMDASolveTExplicit(3)
>>> >>
>>> >> call DMGetLocalVector(daScalars,rhs_loc,ierr);
>>> >> call
>>> DMGlobalToLocalBegin(daScalars,TRHSS_ts,INSERT_VALUES,rhs_loc,ierr);CHKERRQ(ierr)
>>> >> call
>>> DMGlobalToLocalEnd(daScalars,TRHSS_ts,INSERT_VALUES,rhs_loc,ierr);CHKERRQ(ierr)
>>> >
>>> > There is no need to scatter ghost points for rhs_loc.
>>> >
>>> >>
>>> >>
>>> >> call DMRestoreLocalVector(daScalars,t_loc,ierr);CHKERRQ(ierr)
>>> >> call
>>> DMRestoreLocalVector(daScalars,rhs_loc,ierr);CHKERRQ(ierr)
>>> >>
>>> >> end subroutine FormRHSFunction
>>> >>
>>> >> Where DMDASolveTExplicit(3) is the old function to calculate time
>>> integration with runge kutta, modified to only generate the f(t,u) which in
>>> this case is rhs_loc,
>>> >>
>>> >> I still have several doubts:
>>> >>
>>> >> • Will this explicit RHS calculation work with TS routines? my
>>> time integration is explicit so far and it would involve a fair deal of
>>> extra work to make it implicit.
>>> > For explicit time integration, one needs to provide only RHSFunction.
>>> >
>>> >> • This 't' parameter in the RHSFunction is controlled by PETSC?
>>> i am not providing anything for it directly, where is it coming from?
>>> >
>>> > It is controlled by PETSc. If your problem is autonomous (the RHS does
>>> not depend on t), it can be simply ignored.
>>> >
>>> >> • Do i need to provide a Jacobian or the TS routine will try to
>>> guess one? This is related to the first point where my time scheme being
>>> explicit does not use a jacobian.
>>> > For explicit time integration, Jacobian is not needed.
>>> >
>>> > Hong
>>> >
>>> >>
>>> >> Thanks, any help is appreciated, if you see any errors or need more
>>> information please let me know,
>>> >>
>>> >> Regards,
>>> >>
>>> >> Manuel
>>> >>
>>> >> On Wed, Jun 26, 2019 at 9:54 PM Zhang, Hong <hongzhang at anl.gov>
>>> wrote:
>>> >>
>>> >>
>>> >>> On Jun 26, 2019, at 4:17 PM, Manuel Valera via petsc-users <
>>> petsc-users at mcs.anl.gov> wrote:
>>> >>>
>>> >>> Hi PETSc,
>>> >>>
>>> >>> I am trying to implement the Time stepping routines in my model, i
>>> have a working runge-kutta time scheme that goes to the following steps:
>>> >>>
>>> >>> • Interpolate u,v,w to the time advancing variable position.
>>> >>> • Calculate nonlinear coefficients and advect velocities with a
>>> forward-backward shock capturing scheme.
>>> >>> • Calculate the variable laplacian
>>> >>> • Sum terms to create RHS (nonlinear advected velocities +
>>> laplacian)
>>> >>> • Finally, the runge kutta integration is done in a typical way
>>> that looks like:
>>> >>> newtemp(t+1) = prevtemp(t) + dt*RHS
>>> >>>
>>> >>>
>>> >>> So my questions are:
>>> >>> • I think my problem is nonlinear, but is being made linearized
>>> by the advecting scheme, is this correct? this is to know if i should use
>>> the linear or nonlinear routines for TS.
>>> >> TSComputeRHSFunctionLinear is just a convenience function for linear
>>> ODEs in the form udot=Au. Using it won’t buy you much. So for TS starters,
>>> it is fine to assume your problem is nonlinear and think of the form
>>> udot=f(t,u) where f is the RHS function.
>>> >>> • How do i know what are the appropriate routines i should be
>>> using here? so far i think i should use the following:
>>> >>> call TSCreate(PETSC_COMM_WORLD,ts,ierr)
>>> >>> call TSSetProblemType(ts,TS_LINEAR,ierr)
>>> >>> call TSSetTimeStep(ts,dt,ierr)
>>> >>>
>>> >>> call TSSetFromOptions(ts,ierr)
>>> >>>
>>> >>> call TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,NULL,ierr)
>>> >>> call TSSolve(ts,loctemperature,ierr)
>>> >>>
>>> >>> Should i use call TSSetRHSJacobian for the temperature jacobian in
>>> this case?
>>> >>
>>> >> I would suggest to write your own RHSFunction and set it to TS with
>>> TSSetRHSFunction().
>>> >>
>>> >>>
>>> >>>
>>> >>> I am using
>>> https://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/ex4.c.html
>>> as a general guide, is there a more appropriate example?
>>> >>
>>> >> ex4 is a good example. In addition, ex9 uses finite volume method
>>> with slope limiters to solve a variety of problems such as advection
>>> equation, burgers equation and shallow water equation. It might be an
>>> overkill, but it seems to be close to the problem you are trying to solve.
>>> Note that you might want to use the SSP methods (-ts_type ssp) instead of
>>> the classic Runge-Kutta methods for problems with shocks.
>>> >>
>>> >> Hong (Mr.)
>>> >>
>>> >>>
>>> >>> The dt value and total timesteps are controlled by the model,
>>> >>>
>>> >>> Thanks for your help,
>>> >>>
>>> >>>
>>> >>
>>> >
>>>
>>>
>
> --
> 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/
> <http://www.cse.buffalo.edu/~knepley/>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190709/f01324ef/attachment.html>
More information about the petsc-users
mailing list