In `heat.c`, can we replace the SNES solver with KSP solver?

[mapengfei-nwpu]

I have read chapter 5 and I noticed that TS object use the SNES object by default. SNES is a nonlinear solver which calls KSP solver many times in every time step. Is it more expensive than using KSP solver directly?

[mapengfei-nwpu]

In pattern.c, it is obviously that the reaction-diffusion equations can be decoupled into two scalar equations if we take the nonlinear term explicitly. I have tried with PETSc that solving the decoupled equations separately is much faster than solving the whole equations at one time.

[bueler]

Good observations.

Regarding the first point, if you use -snes_type ksponly then the SNES object just calls the KSP once, which is often the right choice if the TS is solving a linear system at each time step. In that case using SNES is not more expensive than KSP. (Using the SNES type newtonls instead might call KSP a second time to clean up solver error, depending on tolerances.) Consistently using SNES means that one can focus on writing a single must-be-correct residual evaluation, and then SNES forms the linear systems without the user needing to write code for it. This is why the examples in the book use SNES.

For your second point, I think what you are describing is exactly what the default ARKIMEX TS type is doing. It treats the nonlinear right-hand-side explicitly and then has two diagonal blocks (decoupled) in the linear system which it solves for the "IFunction" part. So I am curious what you did for "I have tried with PETSc that solving the decoupled equations separately is much faster than solving the whole equations at one time." Can you show the two runs you compared?

[mapengfei-nwpu]

Thank you! I think I didn't describe the second point clearly.

As you write in the book, the nonlinear term is calculated explicitly. Then, what we calculate is a vector-valued heat equation(or a vector-valued Poisson equation in every time step). I think solving a vector-valued Poisson equation is equivalent to solving two scalar-valued Poisson equation.

I modified your program poisson.c to solve a vector-valued function. In function DMDACreate2d, I set the degree of freedom as 2. Finally, the time spending is about 4 times of what scalar-valued Poisson equation spent. Here is my program poisson_vector.txt.

I am very grateful for your respond.

[bueler]

As you write in the book, the nonlinear term is calculated explicitly. Then, what we calculate is a vector-valued heat equation(or a vector-valued Poisson equation in every time step). I think solving a vector-valued Poisson equation is equivalent to solving two scalar-valued Poisson equation.

I think this is a very accurate description of the pattern.c example when the TS type is an IMEX scheme, for example the ARKIMEX default choice. Note that pattern.c does set the degrees of freedom to 2; see line 84.

I modified your program poisson.c to solve a vector-valued function. In function DMDACreate2d, I set the degree of freedom as 2. Finally, the time spending is about 4 times of what scalar-valued Poisson equation spent. Here is my program poisson_vector.txt.

This makes sense for a code which solves two decoupled Poisson equations. But I thought you wanted to solve two coupled diffusion-reaction problems in the time-dependent case? That is what pattern.c is already doing, by the method you want, I think.

[mapengfei-nwpu]

I didn't finished my comment.

Just for Poisson equaiton.

Assume that solving one scalar-valued Poisson equation consumes T seconds. Then, It takes 2T seconds to solve two scalar-valued Poisson equations and 4T seconds to solve one vector-valued Poisson equation. It is tested with poisson_vector and poisson.c. So, I think it is more efficient to solve two scalar-valued equation. This is the core of my opinion.

Can I make statements like this: In default ARKIMEX TS module, the efficency problem that I mentioned have been solved.

[bueler]

In default ARKIMEX TS module, the efficency problem that I mentioned have been solved.

Yes, I think it is solved. That is, for an IMEX scheme on the pattern.c problem the solution of the two scalar poisson problems (at each time step) takes only about twice the time of one scalar Poisson problem. This might be testable by comparing timing for the SNES solver stage of ARKIMEX on the pattern.c problem (dof=2) versus a modified scalar heat equation with some simple nonlinear zeroth-order term (dof=1), again using ARKIMEX. If you do this test then I would be curious to see how it comes out.

[mapengfei-nwpu]

I followed your instruction and did the both tests and get the some results.

for heat.c, I used -ts_monitor -ts_dt 0.005 -ts_adapt_type none -snes_type ksponly , the time consuming: 4.485 for patter.c, I used -ts_monitor -ts_dt 10 -ts_adapt_type none -snes_type ksponly , the time consuming: 13.593 Here, the mesh size is 300x300 and every program solves 20 steps. the time consuming of heat.c is about 3 times of pattern.c.

I quickly modified the pattern.c, and tested it with -ts_monitor -ts_dt 20 -ts_adapt_type none -snes_type ksponly, the time consuming: 3.076. It means original pattern.c(dof=2) takes 4.5 times of the modified pattern.c(dof=1).

If my results is not consistent with yours, I think I should recompile my petsc!