<div dir="ltr"><div dir="ltr">On Tue, Jul 9, 2019 at 2:39 PM Manuel Valera via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr">On Tue, Jul 9, 2019 at 11:27 AM Smith, Barry F. <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">> On Jul 8, 2019, at 6:53 PM, Manuel Valera via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br>
> <br>
> Hi Zhang,<br>
> <br>
> Thanks to your help i have implemented the TS routine for my temperature DMDA array, and it works correctly including in MPI.<br>
> <br>
> 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.<br>
<br>
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. <br></blockquote><div><br></div><div>For some reason this is not possible, if you could provide any ideas let me know, the situation is this:</div><div><br></div><div>This block is now inside the main loop:</div><div> call TSCreate(PETSC_COMM_WORLD,ts,ierr)<br> call TSSetProblemType(ts,TS_NONLINEAR,ierr)<br> call TSSetSolution(ts,LocTemperature,ierr)<br> call TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,PETSC_NULL_VEC,ierr)<br> call TSSetType(ts,TSSSP,ierr)<br> call TSSetTimeStep(ts,dt,ierr)<br> call TSSetDM(ts,daScalars,ierr)<br> call TSSetMaxTime(ts,dt,ierr)<br> call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP,ierr)<br> call TSSetFromOptions(ts,ierr)<br> call TSSetUp(ts,ierr)<br> call TSSolve(ts,LocTemperature,ierr)<br> call TSDestroy(ts,ierr)<br></div><div><br></div><div>Where i create, set everything up and solve the TS problem for a single time step.</div></div></div></blockquote><div><br></div><div>I admit that I have not completely followed the discussion. I will try and ask a few questions so that I can understand.</div><div><br></div><div>First, you can use TSStep(), instead of TSSolve(), in order to run a single time step.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><div> 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</div></div></div></blockquote><div><br></div><div>Is it easy to say what you need to do between steps?</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><div> 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.</div><div><br></div><div>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. </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<br>
> I hope the overhead time is not to bad when i scale to bigger problems.<br>
<br>
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.<br></blockquote><div><br></div><div>Agreed and i would like to know how to reuse them but the behavior seems to be against what i want to do.</div><div> </div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
> <br>
> 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?<br>
<br>
TS doesn't care how you access the vector entries in your function evaluations. All TS every sees is a Vec. <br>
<br>
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.<br>
<br>
But if you mean 3 physical dimensions plus 3 dof then everything is fine just use the 3d DMDA.<br></blockquote><div><br></div><div>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) ? </div></div></div></blockquote><div><br></div><div>There should be no problem with multiple DOF per point in your domain. You just take this into account</div><div>in your residual (and Jacobian) function. There are many examples in TS with multiple DOFs. Maybe there</div><div>is something I am missing here.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> <br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><div>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 ? </div><div><br></div><div>Thanks, i think i am a bit at a loss here,</div><div><br></div><div>Regards,</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<br>
Barry<br>
<br>
<br>
<br>
> <br>
> Thanks, <br>
> <br>
> On Thu, Jul 4, 2019 at 9:18 PM Zhang, Hong <<a href="mailto:hongzhang@anl.gov" target="_blank">hongzhang@anl.gov</a>> wrote:<br>
> <br>
> <br>
>> On Jul 3, 2019, at 3:10 PM, Manuel Valera <<a href="mailto:mvalera-w@sdsu.edu" target="_blank">mvalera-w@sdsu.edu</a>> wrote:<br>
>> <br>
>> Thanks Zhang for your answer,<br>
>> <br>
>> 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,<br>
>> <br>
>> My TS code so far looks like this:<br>
>> <br>
>> (...initialization...)<br>
>> <br>
>> call TSCreate(PETSC_COMM_WORLD,ts,ierr)<br>
>> call TSSetProblemType(ts,TS_NONLINEAR,ierr)<br>
>> call TSSetSolution(ts,gTemperature,ierr)<br>
>> call TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,PETSC_NULL_VEC,ierr)<br>
> <br>
> The second last argument should be the user context. If you pass NULL, the ctx variable in your FormRHSFunction will be NULL. <br>
> <br>
>> call TSSetType(ts,TSSSP,ierr)<br>
>> call TSSetTimeStep(ts,dt,ierr)<br>
>> call TSSetDM(ts,daScalars,ierr)<br>
>> call TSSetMaxTime(ts,TotalTime,ierr)<br>
>> call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)<br>
>> call TSSetFromOptions(ts,ierr)<br>
>> call TSSetUp(ts,ierr)<br>
>> <br>
>> (... then inside main calculations loop i have...)<br>
>> <br>
>> call TSSolve(ts,gTemperature,ierr)<br>
>> <br>
>> (...and my RHSFunction looks like this...)<br>
>> <br>
>> subroutine FormRHSFunction(ts,t,t_loc,rhs_loc,ctx,ierr)<br>
>> <br>
>> use petscvec<br>
>> use petscts<br>
>> use petscdmda<br>
>> use petscdm<br>
>> <br>
>> USE utils<br>
>> USE dmdaobjs<br>
>> USE probsize<br>
>> USE modelparams<br>
>> <br>
>> implicit none<br>
>> <br>
>> TS :: ts<br>
>> PetscReal :: t<br>
>> Vec :: ctx,t_loc,rhs_loc<br>
>> PetscErrorCode :: ierr<br>
>> <br>
>> <br>
>> <br>
>> call TSGetDM(ts,daScalars,ierr)<br>
>> <br>
>> call DMGetLocalVector(daScalars,t_loc,ierr);<br>
> <br>
> t_loc is the input, it should be not modified.<br>
> <br>
>> call DMGlobalToLocalBegin(daScalars,gTemperature,INSERT_VALUES,t_loc,ierr);CHKERRQ(ierr)<br>
>> call DMGlobalToLocalEnd(daScalars,gTemperature,INSERT_VALUES,t_loc,ierr);CHKERRQ(ierr)<br>
> <br>
> FormRHSFunction is supposed to implement rhs_loc = f(t,t_loc). So you want to scatter ghost points for t_loc, not gTemperature.<br>
> <br>
>> call DMDASolveTExplicit(3)<br>
>> <br>
>> call DMGetLocalVector(daScalars,rhs_loc,ierr);<br>
>> call DMGlobalToLocalBegin(daScalars,TRHSS_ts,INSERT_VALUES,rhs_loc,ierr);CHKERRQ(ierr)<br>
>> call DMGlobalToLocalEnd(daScalars,TRHSS_ts,INSERT_VALUES,rhs_loc,ierr);CHKERRQ(ierr)<br>
> <br>
> There is no need to scatter ghost points for rhs_loc.<br>
> <br>
>> <br>
>> <br>
>> call DMRestoreLocalVector(daScalars,t_loc,ierr);CHKERRQ(ierr)<br>
>> call DMRestoreLocalVector(daScalars,rhs_loc,ierr);CHKERRQ(ierr)<br>
>> <br>
>> end subroutine FormRHSFunction<br>
>> <br>
>> 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,<br>
>> <br>
>> I still have several doubts:<br>
>> <br>
>> • 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. <br>
> For explicit time integration, one needs to provide only RHSFunction.<br>
> <br>
>> • This 't' parameter in the RHSFunction is controlled by PETSC? i am not providing anything for it directly, where is it coming from?<br>
> <br>
> It is controlled by PETSc. If your problem is autonomous (the RHS does not depend on t), it can be simply ignored.<br>
> <br>
>> • 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.<br>
> For explicit time integration, Jacobian is not needed.<br>
> <br>
> Hong<br>
> <br>
>> <br>
>> Thanks, any help is appreciated, if you see any errors or need more information please let me know,<br>
>> <br>
>> Regards,<br>
>> <br>
>> Manuel <br>
>> <br>
>> On Wed, Jun 26, 2019 at 9:54 PM Zhang, Hong <<a href="mailto:hongzhang@anl.gov" target="_blank">hongzhang@anl.gov</a>> wrote:<br>
>> <br>
>> <br>
>>> On Jun 26, 2019, at 4:17 PM, Manuel Valera via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br>
>>> <br>
>>> Hi PETSc,<br>
>>> <br>
>>> 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:<br>
>>> <br>
>>> • Interpolate u,v,w to the time advancing variable position.<br>
>>> • Calculate nonlinear coefficients and advect velocities with a forward-backward shock capturing scheme.<br>
>>> • Calculate the variable laplacian <br>
>>> • Sum terms to create RHS (nonlinear advected velocities + laplacian)<br>
>>> • Finally, the runge kutta integration is done in a typical way that looks like:<br>
>>> newtemp(t+1) = prevtemp(t) + dt*RHS<br>
>>> <br>
>>> <br>
>>> So my questions are:<br>
>>> • 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.<br>
>> 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. <br>
>>> • How do i know what are the appropriate routines i should be using here? so far i think i should use the following:<br>
>>> call TSCreate(PETSC_COMM_WORLD,ts,ierr)<br>
>>> call TSSetProblemType(ts,TS_LINEAR,ierr)<br>
>>> call TSSetTimeStep(ts,dt,ierr)<br>
>>> <br>
>>> call TSSetFromOptions(ts,ierr)<br>
>>> <br>
>>> call TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,NULL,ierr)<br>
>>> call TSSolve(ts,loctemperature,ierr)<br>
>>> <br>
>>> Should i use call TSSetRHSJacobian for the temperature jacobian in this case?<br>
>> <br>
>> I would suggest to write your own RHSFunction and set it to TS with TSSetRHSFunction().<br>
>> <br>
>>> <br>
>>> <br>
>>> I am using <a href="https://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/ex4.c.html" rel="noreferrer" target="_blank">https://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/ex4.c.html</a> as a general guide, is there a more appropriate example?<br>
>> <br>
>> 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.<br>
>> <br>
>> Hong (Mr.)<br>
>> <br>
>>> <br>
>>> The dt value and total timesteps are controlled by the model,<br>
>>> <br>
>>> Thanks for your help,<br>
>>> <br>
>>> <br>
>> <br>
> <br>
<br>
</blockquote></div></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>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</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>