Class DenseOutput to get dense output of solution
#464
Conversation
…ing some documentation.
There was a problem hiding this comment.
Seems a bit more complicated than needed in places. One of my favourite design principles is KISS, which always reminds me that we are all just a bunch of idiots. ;)
To me, KISS applies mainly in refactoring. After I "finish" something to the point that it does what I want, I go back and check how much more it is doing that I need and how much less generic it is than possible. This is best done by unit testing which tests only a single feature rather then the whole thing. This often exposes what input is really needed for a function. If you think it's needed, you can add an integration test afterwards to test the whole thing is working.
Everybody has their own style of coding. I learnt that writing the unit tests at the same time as I write the functions is what works best for me. As Robert's book recommendation The Pragmatic Programmer points out, the test is the first client of your code. Meaning that as you switch between writing code and test, you switch between being developer and user. Crucially, you should put yourself in the role of a generic user. Don't develop the test for exactly what you need, do it as generically as possible.
| @@ -6,19 +6,25 @@ | |||
|
|
|||
| class LogSolution(Hooks): | |||
| """ | |||
| Store the solution at the end of each step as "u". | |||
| Store the solution at the end of each step as "u". It is also possible to get the solution | |||
There was a problem hiding this comment.
Please put this in a separate hook. For large simulations, this may be a huge amount of data and possibly a huge amount of communication. This should not be the default behaviour.
| uTmp = L.u if not L.sweep.coll.left_is_node else L.u[1:] | ||
|
|
||
| # Determine the shape of arrays to handle None values | ||
| shape = max((me.shape for me in uTmp if me is not None), default=(0,)) |
There was a problem hiding this comment.
This is fishy to me. Why would the entries have different shapes? Why not simply use the shape of L.u[0], which all ranks have?
| recvBuf = np.empty(size * len(uFlat), dtype='d') | ||
|
|
||
| # Use Allgather to collect data from all processes | ||
| comm.Allgather(uFlat, recvBuf) |
There was a problem hiding this comment.
Do we need to gather on all ranks, actually? Why not gather it on rank 0 only?
| uFlat = np.concatenate(u) | ||
|
|
||
| # Prepare the buffer to receive data from all processes | ||
| recvBuf = np.empty(size * len(uFlat), dtype='d') |
There was a problem hiding this comment.
what about recvBuf = np.empty_like(uFlat)? If you specify d as datatype, this cannot be used for complex data.
| comm.Allgather(uFlat, recvBuf) | ||
|
|
||
| # Reshape and combine data from all processes | ||
| recvBuf = recvBuf.reshape(size, -1, shape[0]) |
There was a problem hiding this comment.
Do we need to communicate in flattened form? I think mpi4py can handle any numpy array as buffer, no?
| """Initialization routine""" | ||
|
|
||
| # Extract the stage times and values | ||
| self.nodes = [np.array(entry[1]) for entry in nodes] |
There was a problem hiding this comment.
please use numpy arrays here. I guess your 40 days simulation will generate 1e9123832 data points. The difference between numpy arrays and lists will likely be around 14 ages of the universe. (Rough estimate ;))
| """ | ||
|
|
||
| for index, tau in enumerate(self.nodes): | ||
| if tau[0] <= t <= tau[-1]: |
There was a problem hiding this comment.
You should use masks here for faster computation.
| uValuesSub = [subElement[subIndex] for subElement in uValuesComponents] | ||
|
|
||
| # Create an LagrangeApproximation object for this dimension | ||
| sol = LagrangeApproximation(points=nodes, fValues=uValuesSub) |
There was a problem hiding this comment.
Do you need this for all components? The interpolation matrix should be the same, no? Why not generate one interpolation matrix and multiply all component values?
There was a problem hiding this comment.
Why do you need to worry about the components anyways? Multiplying the 1xM interpolation matrix with an MxXxYxZ matrix return an XxYxZ array, so you can interpolate everything in one operation with no regard as to what kind of mesh this is, no?
| @pytest.mark.base | ||
| @pytest.mark.parametrize("quad_type", ['RADAU-RIGHT', 'LOBATTO']) | ||
| @pytest.mark.parametrize("problemType", ['ODE', 'DAE']) | ||
| def test_interpolate(quad_type, problemType): |
There was a problem hiding this comment.
I would prefer you not generate the data using a simulation. Ideally, tests will target as small a unit of code as possible. So the tests should fail only if what they test breaks not if somebody changed some parameter somewhere else or whatever. Also, it's faster to just generate a dictionary than to run a simulation.
|
|
||
| import numpy as np | ||
|
|
||
| nodesRandom = { |
There was a problem hiding this comment.
Well, it's not really random... This is overselling a bit :D
Totally agree with you ... and I would suggest something similar to your messages @brownbaerchen 😉 |
As already mentioned in the pySDC v6 discussion it would be nice to get the solution also at points not computed in the simulation by doing interpolation. At least I need that for computing a reference solution of the Van der Pol DAE problem that has a quite steep slope requiring a very small time step size$\Delta t\leq 10^{-9}$ to have a good resolution of this. However the computation takes forever.. Thus, a computation of adaptivity to speedup the computation might make a sense here!
To realize this I implemented hooks to log the solution at all nodes together with nodes in
LogSolutionand provided a test for that. After simulation you obtainand then the class
DenseOutputis used:Evaluating at requested times
t_evalyou getThe class
DenseOutputrealizes the handling of calling the numerical solution and interpolation at requested time points. For interpolationqmatpackage is used.I was not sure where to store the
DenseOutputclass inpySDCso I decided to first store it in theDAEproject together with the test.I assume it is a nice-to-have here also in$1$
pySDCeven routines likescipy.integrate.solve_ivpalso have that feature. For my issue: Since the Van der Pol DAE is a DAE of indexscipy.integrate.solve_ivpcan still be used for solving. However, due to the steep slope the routine fails at some point and provided the error message thatRequired step size is less than spacing between numberswhere for this problem solving this very stiff example seem to not suitable anymore. Nice to know: I discovered that someone is working on a routine of scipy to solve DAEs. I was very happy about that and played around with that. So I put Van der Pol into it and look! I got the same error message.. so sad.. Anyways. Using theDenseOutputthis should speedup the things.What do you think about that? As usual: Improvements, corrections, comments, etc. are really appreciated! Many thanks in advance! :)