Added solve_system to all DAE problem classes

[lisawim]

Currently, the fully_implicit_DAE sweeper solves all DAE problems in the same way, i.e., using scipy.optimize.root(f, x0, method='hybr'). Since it could be useful to choose the solver problem dependent, a solve_system method is added to all DAE classes as this is already the case for all implemented ODE/PDE classes in pySDC. The implicit system is defined in the sweeper anyway, because this nonlinear system does not change (since it is the collocation problem to be solved).

So no magic here in this PR, I only moved these two line

opt = root(impl_sys, u0, method='hybr', tol=self.newton_tol)
me[:] = opt.x

Now, the problem classes still uses the same solver, but for the future, also different solver could be implemented and no one has to take care to change the sweeper itself every time.
to the problem classes.

Now, also WorkCounters are used to count the work which is also useful to get some statistics about the calls of the right-hand side and the calls of the right-hand side during the solve (because for hybr no callback is available and also no iterations needed for the solve can be found..).

[pancetta]

Very nice! To me there seems to be quite a bit of boilerplate code (since the solve_system routine is there for many different problem classes. Would it make sense to have this as a generic/default solver routine somewhere? If there is even more boilerplate code like this (didn't check), you could inherit from one of the problems classes when designing the other ones?

[lisawim]

Of course it does make sense to have a class with the solve_system method and all problem classes could inherit from! I was thinking of the pySDC.projects.DAE.misc.ProblemDAE class where the method is added, because all DAE classes inherits from it anyway. Would this be okay? Since another solver should be implemented, solve_system can then be overwritten in the corresponding class.

[pancetta]

Sounds good! You can add more generic functions there, I assume?

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (73cb8e3) 73.62% compared to head (ff5805f) 73.63%.

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #371   +/-   ##
=======================================
  Coverage   73.62%   73.63%           
=======================================
  Files         270      270           
  Lines       22383    22380    -3     
=======================================
  Hits        16480    16480           
+ Misses       5903     5900    -3     

Files	Coverage Δ
pySDC/projects/DAE/misc/ProblemDAE.py	100.00% <100.00%> (ø)
pySDC/projects/DAE/problems/simple_DAE.py	97.87% <100.00%> (+0.14%)	⬆️
pySDC/projects/DAE/problems/synchronous_machine.py	98.63% <100.00%> (+0.01%)	⬆️
...ySDC/projects/DAE/problems/transistor_amplifier.py	96.29% <100.00%> (+0.14%)	⬆️
pySDC/projects/DAE/sweepers/fully_implicit_DAE.py	97.36% <100.00%> (+2.68%)	⬆️

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[lisawim]

Yes, this should be possible.

[pancetta]

Looks like at least this one could go to the generic parent class, too?

[lisawim]

True! I added this for all problem classes. I'll do that.

[lisawim]

Looks better and more convenient (like before), thanks a lot!

[brownbaerchen]

To me it seems like you could indeed save a lot of lines by inheriting from generic_implicit in the fully_implicit_DAE sweeper. You can also consider doing multiple inheritance if you plan to do IMEX at some point.
I have not checked very thoroughly, but I believe the integrate, compute_residual, predict and compute_end_point functions are all equivalent to the ones in generic_implicit.
If you move the implSystem function to the problem class, would it be possible to just use the generic_implicit sweeper? I think there was some issue with this. Also, u_approx is not usually passed to solve_system, I think..? It probably doesn't work, but if it did, that would be nice.

[brownbaerchen]

decides -> decide

[lisawim]

To me it seems like you could indeed save a lot of lines by inheriting from generic_implicit in the fully_implicit_DAE sweeper. You can also consider doing multiple inheritance if you plan to do IMEX at some point. I have not checked very thoroughly, but I believe the integrate, compute_residual, predict and compute_end_point functions are all equivalent to the ones in generic_implicit. If you move the implSystem function to the problem class, would it be possible to just use the generic_implicit sweeper? I think there was some issue with this. Also, u_approx is not usually passed to solve_system, I think..? It probably doesn't work, but if it did, that would be nice.

Yes, IMEX is indeed planned (and already done but not pushed yet), and I agree that multiple inheritance can be done here. At least integrate, compute_end_point and predict could be used from generic_implicit. For compute_residual I have to think about that. eval_f is used here anyway, but in case of DAEs we also have the du argument here. This issue could probably be solved by setting du as a arbitrary argument (i.e., by setting du=None somewhere)..?
Moving the implSystem function to the problem class - I'll check if this would work.

[lisawim]

Hm, I also got a message that you @brownbaerchen wrote somewhere (I cannot find it) that the scipy.optimize.root routine has a callback where I can call work_counters['newton']. This is true, but for method='hybr' this callback function is not available.

[brownbaerchen]

Hm, I also got a message that you @brownbaerchen wrote somewhere (I cannot find it) that the scipy.optimize.root routine has a callback where I can call work_counters['newton']. This is true, but for method='hybr' this callback function is not available.

Yes, I saw what you wrote after I made the comment, so I marked it as resolved.

[brownbaerchen]

To me it seems like you could indeed save a lot of lines by inheriting from generic_implicit in the fully_implicit_DAE sweeper. You can also consider doing multiple inheritance if you plan to do IMEX at some point. I have not checked very thoroughly, but I believe the integrate, compute_residual, predict and compute_end_point functions are all equivalent to the ones in generic_implicit. If you move the implSystem function to the problem class, would it be possible to just use the generic_implicit sweeper? I think there was some issue with this. Also, u_approx is not usually passed to solve_system, I think..? It probably doesn't work, but if it did, that would be nice.

Yes, IMEX is indeed planned (and already done but not pushed yet), and I agree that multiple inheritance can be done here. At least integrate, compute_end_point and predict could be used from generic_implicit. For compute_residual I have to think about that. eval_f is used here anyway, but in case of DAEs we also have the du argument here. This issue could probably be solved by setting du as a arbitrary argument (i.e., by setting du=None somewhere)..? Moving the implSystem function to the problem class - I'll check if this would work.

In this case, I vote to just overload the function and copy paste a few lines rather than adding default arguments to the other stuff. Multiple inheritance is tricky business in my opinion. For me, it requires a bit of trial and error as well as some praying. But it is kind of working here and I think it's quite elegant once it works.

[lisawim]

Indeed, I could removed lines by removing integrate and compute_end_point in the sweeper. For predict, compute_residual I didn't see yet how I can use the one from the parent class. For predict using spread eval_f is used here to spread the values. For compute_residual also eval_f is used to define the residual.

[brownbaerchen]

I see. In predict, I think the only thing that's different is spread. So you could call the parent's spread method if the initial guess is not spread and otherwise do what you have here. But that's optional :D

As for the residual, I am a bit confused that the residual is based on the magnitude of the right hand side evaluations rather than the integral. But I have no idea about the DAE implementation. Feel free to ignore this. You don't have to abstract as much as you can. Using code that is already there is better than repeating code, but bending over backwards to somehow make the code that is already there do what you want in a different place is also no good.

[lisawim]

I tried that but in the predict method from sweeper it is L.f[0] = P.eval_f(L.u[0], L.time) used to compute the function value at initial time, and here I also have the struggle by using eval_f which needs du to be computed..

Regarding the residual: At first I was also a bit confused because using the right-hand side for computing the residual was a bit cumbersome. Let me explain you shortly (as this makes sense for me at least): How do you compute the residual for a system of equations? In case of having $Au=b$, the residual can be computed by $r=Au-b$. In case of nonlinearity, you can define a system like $F(u)=0$, then $r=F(u)$ is the residual. When you think of that a DAE is given in the general nonlinear formulation $$F\left(\frac{d u}{dt}, u, t\right) = 0,$$ then you can also compute the residual by evaluating the right-hand side using your current numerical solutions L.u and L.f here. So, for the current numerical solution it is checked how well it does satisfy the DAE formulation or not. Does that answer your question?

But now I'm thinking whether we also can compute the residual by using the integral formulation inside the DAE formulation.. Hm..

I know that repeating code is bad practice. @pancetta have already taught me the DRY principle some time, but I don't always see where I can apply it. In this case, I'm very thankful for such comments!

[lisawim]

Ah, I see. For $F\left(u, \frac{du}{dt}, t\right)=0$, you have the (node-wise) collocation formulation $$F\left(u_0 + \sum_{j=1}^M q_{m, j} U_j, U_m, \tau_m\right)=0, \qquad m=1,..,M$$ for $u_m\approx u(\tau_m)$, $U_m \approx u'(\tau_m)$ and collocation nodes $t_0 \leq \tau_1 <..<\tau_M \leq t_1$.
This would also be the formula to compute the (node-wise) residual in case of DAEs. The integral is already computed in line 148-149 in fully_implicit_DAE, so you can use L.u[m + 1] here, insert into eval_f and this is the residual at $\tau_m$. Does that make sense?

[brownbaerchen]

Honestly, don't worry about abstraction. I was missing some small details that do make a difference and I think there is little to gain from further abstraction. Sorry about that :D

[pancetta]

Can we merge this? Issue lisawim#8 seems to be a topic for another PR anyway?

[lisawim]

Yes, can be merged.