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 <[email protected]> wrote: > On Tue, Jul 9, 2019 at 2:39 PM Manuel Valera via petsc-users < > [email protected]> wrote: > >> On Tue, Jul 9, 2019 at 11:27 AM Smith, Barry F. <[email protected]> >> wrote: >> >>> > On Jul 8, 2019, at 6:53 PM, Manuel Valera via petsc-users < >>> [email protected]> 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 <[email protected]> wrote: >>> > >>> > >>> >> On Jul 3, 2019, at 3:10 PM, Manuel Valera <[email protected]> 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 <[email protected]> >>> wrote: >>> >> >>> >> >>> >>> On Jun 26, 2019, at 4:17 PM, Manuel Valera via petsc-users < >>> [email protected]> 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/> >
