# PETSc example (heat-TS-2d-PETSc.edp) incorrectly imlpemented? **URL:** \ **Category:** General Discussion\ **Created:** [January 23, 2025, 7:32am UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704 "2025-01-23T07:32:51Z")\ **Posts on this page:** 13\ **Page:** 1
**Author:** ![usiu5555](https://avatars.discourse-cdn.com/v4/letter/u/e19adc/32.png) [@usiu5555](https://community.freefem.org/u/usiu5555)\ **Post date:** [January 23, 2025, 7:32am UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/1 "2025-01-23T07:32:51Z")
I am working with the [example](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/heat-TS-2d-PETSc.edp). I just tried to change the boundary condition (line 24, on(1, u=0)) to a different value, and the solution does not change. It seems like something is missing, as regardless of the value substituted with 0, the solution always is the same as for u=0. Initially, I wanted to include Neumann BC (int2d(Th,1)(v\*q)), but that did not work. Any advice on how can I apply different boundary conditions while using PETSc? ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [January 23, 2025, 9:11am UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/2 "2025-01-23T09:11:09Z")
You need to adjust the initial guess. ---
**Author:** ![usiu5555](https://avatars.discourse-cdn.com/v4/letter/u/e19adc/32.png) [@usiu5555](https://community.freefem.org/u/usiu5555)\ **Post date:** [January 24, 2025, 3:29am UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/4 "2025-01-24T03:29:06Z")
I don’t understand why the initial condition should affect the boundary conditions. For example, it should work if I have an initial condition equal to 1 everywhere and then at t=0 demand value 10 at a given boundary. Indeed, this works when the value is set to 0, but only to 0. This should work for any initial condition with Neumann’s boundary condition as well. Could you please elaborate? What I found is that the example heat-TS-RHS-2d-PETSc.edp is probably the key to the problem. Setting ```auto w[] = rhs; ``` inside funcRHS seems to impose the boundary conditions, but definitely not in the right way. Is that the correct approach to the problem? Even though something is wrong, it seems to be closer to the solution. ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [January 24, 2025, 11:30am UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/5 "2025-01-24T11:30:11Z")
> [@usiu5555](#): > > I don’t understand why the initial condition should affect the boundary conditions I slightly fixed `heat-TS-2d-PETSc.edp`: ```auto diff --git a/examples/hpddm/heat-TS-2d-PETSc.edp b/examples/hpddm/heat-TS-2d-PETSc.edp index 2a31e596..a84b4511 100644 --- a/examples/hpddm/heat-TS-2d-PETSc.edp +++ b/examples/hpddm/heat-TS-2d-PETSc.edp @@ -22,7 +22,7 @@ matrix Loc; // local operator fespace Ph(Th, P0); Ph kappa = x < 0.25 ? 10.0 : 1.0; varf vPb(u, v) = int2d(Th)(-1.0 * kappa * grad(u)' * grad(v)) + on(1, u = 0.0); - Loc = vPb(Wh, Wh, tgv = -2); + Loc = vPb(Wh, Wh, tgv = -10); rhs = vPb(0, Wh, tgv = -2); } ``` Now, if you switch to something else than `0` on the boundary, e.g., `w = (0.5 - x)^2 + (0.5 - y)^2 < 0.2 ? 2.0 : 1.0;`, you’ll should have the proper conditions. ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [January 24, 2025, 2:36pm UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/6 "2025-01-24T14:36:38Z")
Apart from boundary conditions and right-hand side difficulties in this code (`rhs` is nowhere used), it lacks the mass lumped matrix. Here is a correction with an analytic solution as test case. [heat-TS-2d-PETSc.edp](https://community.freefem.org/uploads/short-url/inp3s0aYvw6r2aBtg7K5ezlnppW.edp) (3.5 KB) ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [January 26, 2025, 11:15pm UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/7 "2025-01-26T23:15:30Z")
A better implementation with full mass matrix [heat-TS-2d-PETSc.edp](https://community.freefem.org/uploads/short-url/jZX4hm01xWjPiLK0V14AErolqFw.edp) (2.9 KB) ---
**Author:** ![usiu5555](https://avatars.discourse-cdn.com/v4/letter/u/e19adc/32.png) [@usiu5555](https://community.freefem.org/u/usiu5555)\ **Post date:** [January 29, 2025, 4:01am UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/8 "2025-01-29T04:01:38Z")
Thank you for the update. Can you explain why changing this line was necessary? Also, the boundary condition is imposed on line 23 (one line above the one you changed, the term on(1, u = 0.0)), so changing it through the initial conditions seems off. It also doesn’t allow to set Neumann’s BC. Looking at fb77 reply, it seems like the problem indeed was deeper. ---
**Author:** ![usiu5555](https://avatars.discourse-cdn.com/v4/letter/u/e19adc/32.png) [@usiu5555](https://community.freefem.org/u/usiu5555)\ **Post date:** [January 29, 2025, 4:07am UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/9 "2025-01-29T04:07:51Z")
This looks like what I was looking for, thank you 🙂 I think I will be able to understand most of the changes and Neumann boundary condition also works. If I am correct, the funcRes was the problem, and also tgv value. Could you comment on the parameter tgv and why do you set the specific values (-1 for rhs and -10 for mass matrix)? Also, at my goal I woud like to impose nonlinear boundary conditions (flux at the boundary as a nonlinear function of u). Will I be able to achieve that with the use of TSSolve? ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [January 29, 2025, 11:01am UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/10 "2025-01-29T11:01:08Z")
There is a document describing PETSc/TS at > **[PETSc/TS: A Modern Scalable ODE/DAE Solver Library](https://arxiv.org/abs/1806.01437)** > > High-quality ordinary differential equation (ODE) solver libraries have a long history, going back to the 1970s. Over the past several years we have implemented, on top of the PETSc linear and nonlinear solver package, a new general-purpose,... Looking at the array p. A:11, the line that is relevant for us is line 3 F(t,u,\dot u):=M\dot u-h(t,u), meaning that TS solves the differential equation F(t,u,\dot u)=0. Here M is the mass matrix. This F corresponds to `funcRes`, the argument `in` corresponds to u, and the argument `inT` corresponds to \dot u. The application of Dirichlet BC is a bit subtle. It is of the form u\_j=u\_j^D(t), where u\_j stands for a particular component of the vector u (a component on the Dirichlet boundary), and u\_j^D(t) is the given Dirichlet data. This equation is not a differential equation in time. The equation on u\_j is encoded in the component j of F. Thus the line j of M has to be set to zero, and one has to set -h\_j(t,u):=u\_j-u\_j^D(t) This is what is done with `tgv=-10` for M and `tgv=-1` for h\_j, see the FreeFem documentation p.221/222 > **[FreeFEM-documentation.pdf](https://doc.freefem.org/pdf/FreeFEM-documentation.pdf)** > > 33.85 MB The TS procedure is able to solve nonlinear equations (F is nonlinear). I think that this is done by Newton’s method. This is why we have to provide the Jacobian matrix in `funcJ` (see p. A:4). If you want to solve nonhomogeneous Neumann condition with a nonlinearity, I think this is possible, you have to put it in the variational formulation `vPbrhs` and using `in` in place of u (but don’t forget to do `ChangeNumbering` each time you switch from Freefem vectors to PETSc vectors or the converse). You have also to modify the Jacobian matrix and make it dependent on `in`. About the time integration method, instead of `-ts_type beuler`, (backward Euler) which is 1rst order, I recommand to use `cn` (Crank Nicholson) or `bdf`, which are second order. I experienced that it is interesting to take constant timestep by setting `-ts_adapt_type none` (otherwise, timestep adaptation with `atol` and `rtol` is described p. A:17) Francois. ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [January 29, 2025, 1:33pm UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/11 "2025-01-29T13:33:03Z")
You may prefer a handmade time loop like this [heat-PETSc.edp](https://community.freefem.org/uploads/short-url/iu2ckdcRubpMCMb8MeHqnXEsPgr.edp) (1.7 KB) Then you have to do second-order, Newton, timestep adaptation by hand. But it is much more understandable than with TSSolve. ---
**Author:** ![usiu5555](https://avatars.discourse-cdn.com/v4/letter/u/e19adc/32.png) [@usiu5555](https://community.freefem.org/u/usiu5555)\ **Post date:** [February 3, 2025, 4:19am UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/12 "2025-02-03T04:19:52Z")
Thank you for the detailed answer and the second script, this helps a lot. In the end, I want to switch to CN, but what do you mean by “it is interesting to take constant timestep”? In my experience using timestep adaptation is a great safeguard to ensure appropriate timestep throughout the simulation. I am surprised that the initial timestep is not checked by the TSS and is accepted regardless of the associated error. Indeed I am more familiar with the handmade loops and I might be forced to do my project this way, but I hope I will be able to stick to the professional implementation. ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [February 3, 2025, 8:33am UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/13 "2025-02-03T08:33:24Z")
I don’t know exactly how TS manages the timestep, in particular what it does with the one given by the user via `-ts_dt` (except for `-ts_adapt_type none` in which case it just takes this timestep). I agree with you on the usefullness of timestep adaptation in order to save computational time. But I noticed that taking constant timestep in TS gives a really smaller error than with adapted timestep. Thus I am wondering if there could be a bug in the TS implementation, that would lower the order of accuracy for non constant timestep (otherwise we have to understand how to choose the parameters `atol` or `rtol` so as to save computational time without loosing accuracy) . Using a handmade time loop allows us to be sure of what we are doing! ---
**Author:** ![usiu5555](https://avatars.discourse-cdn.com/v4/letter/u/e19adc/32.png) [@usiu5555](https://community.freefem.org/u/usiu5555)\ **Post date:** [February 5, 2025, 6:44am UTC](https://community.freefem.org/t/petsc-example-heat-ts-2d-petsc-edp-incorrectly-imlpemented/3704/14 "2025-02-05T06:44:04Z")
Interesting. Indeed TS implementation lacks some control and shows undesired behavior (like the error you mentioned and the initial timestep not being tested). I agree with your point, especially that implementing a basic timestep control is not very complicated. However, I prefer to use professionally implemented solutions if possible. That’s why I started using FreeFEM 🙂 Thank you for all the help on this topic, I hope it will be useful for others and maybe your solution can replace, or be used to improve, the original example.