I am working with the example.
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?
You need to adjust the initial guess.
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
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.
I slightly fixed heat-TS-2d-PETSc.edp:
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<real> 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.
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 (3.5 KB)
A better implementation with full mass matrix
heat-TS-2d-PETSc.edp (2.9 KB)
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.
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?
There is a document describing PETSc/TS at
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
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.
You may prefer a handmade time loop like this
heat-PETSc.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.
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.
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!
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.