parallel-in-time / pysdc Goto Github PK

Modify test_imexsweeper so that it tests all implemented types of quadrature nodes.

Right now, the functions to generate the stability function and matrix representation for SDC are copy&pasted at different places (tests, plot_stability). They should all me moved to one place.

Looks like the new collocation classes still have issues:
https://github.com/Parallel-in-Time/pySDC/runs/5582155629?check_suite_focus=true#step:9:407
Any ideas, @tlunet ?

Implement tests for

• Sweeper
• imex_1st_order

and make sure that the derived update matrix gives the same result as performing k many sweeps.

This will allow a rudimentary multi-level hierarchy and some tests for fault tolerance.

It seems that the errors are wrong for all Tend < 0.5. If I do just one step with dt=0.125, the error is very high. The same is true for 2 and 3 steps, while for 4 steps all is good..

When trying to run a pySDC example with
Python 2.7.10 :: Anaconda 2.3.0 (64-bit)
I get the following error message from the future package:

Traceback (most recent call last): File "playground.py", line 10, in <module> from pySDC.sweeper_classes.imex_1st_order import imex_1st_order File "/usr/not-backed-up/pySDC/pySDC/sweeper_classes/imex_1st_order.py", line 3, in <module> from pySDC.Sweeper import sweeper File "/usr/not-backed-up/pySDC/pySDC/Sweeper.py", line 7, in <module> from future.utils import with_metaclass ImportError: No module named future.utils

Integrate the code provided by O. Knoth via email into standard_integrators.

Can we integrate the function check_iteration_estimator into the new convergence_controllers? Maybe I'm wrong, but it looks very similar. What do you think, @brownbaerchen?

For the fast-wave slow-wave problem at least, it seems that there is no difference between setting maxiter=0 and maxiter=1 (for SDC, using num_proc=1) since reported errors are identical in all digits.

Also, the stability plot in both cases is different from that of IMEX Euler: it should be stable up to 2.0 on the slow axis, because dt=1.0 and forward Euler is stable until 2.0/abs(a) (cf. attached plot).

Is this expected behaviour?

NOTE: Compared wrong things here, plot showed standard IMEX, not fast-wave slow-wave. For the latter, there are still differences but they are much smaller.

Even for Gauss-Lobatto, a simple copy from u_M to u_end is not correct, if you are picky. To reproduce a RK/collocation scheme you will have to evaluate the collocation formula in ALL cases, not just for Gauss-Legendre et al.. Fix this e.g. here:
https://github.com/Parallel-in-Time/pySDC/blob/master/pySDC/sweeper_classes/imex_1st_order.py#L146

Does this change anything...?

When trying to run the penningtrap_playground.py I get the following:

Traceback (most recent call last):
File ".../pySDC-2.1/playgrounds/Boris/penningtrap_playground.py", line 89, in
main()
File ".../pySDC-2.1/playgrounds/Boris/penningtrap_playground.py", line 72, in main
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
File ".../pySDC-2.1/pySDC/implementations/controller_classes/allinclusive_classic_nonMPI.py", line 80, in run
self.hooks.pre_run(step=S, level_number=0)
File ".../pySDC-2.1/playgrounds/Boris/penningtrap_HookClass.py", line 46, in pre_run
N = L.prob.nparts
AttributeError: 'penningtrap' object has no attribute 'nparts'

There seems to be an issue with interpolation between quadrature nodes if the difference in the number of nodes between levels is too large.

Running e.g. Step 5 from the tutorial with

sweeper_params['num_nodes'] = [7, 3]

terminates with message

AssertionError: ERROR: error is too high, got 0.000420431707947

In contrast,

sweeper_params['num_nodes'] = [3, 3]

works fine.

To 1D acoustic-advection and Boussinesq, add a DIRK method of order 5 to standard_integrators.

Coefficients can be found here: http://link.springer.com/article/10.1007%2FBF01990538#page-1

Currently trying to run and compare with the Boris playgrounds.
'penningtrap_playground.py' runs fine when invoked from the terminal or IPython console, but 'spiraling_particle_playground.py' throws a value errror when run using any of those methods:

` uinit, stats, a0 = main(dt=1.0, Tend=5000)

File "E:\Documents\OneDrive - University of Leeds\phd\code\pySDC-master\pySDC\playgrounds\Boris\spiraling_particle_playground.py", line 68, in main
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)

File "C:\ProgramData\Anaconda3\lib\site-packages\pySDC\implementations\controller_classes\allinclusive_classic_nonMPI.py", line 105, in run
self.MS[p] = self.pfasst(self.MS[p], len(active_slots))

File "C:\ProgramData\Anaconda3\lib\site-packages\pySDC\implementations\controller_classes\allinclusive_classic_nonMPI.py", line 310, in pfasst
S.levels[0].sweep.compute_residual()

File "C:\ProgramData\Anaconda3\lib\site-packages\pySDC\core\Sweeper.py", line 169, in compute_residual
res = self.integrate()

File "C:\ProgramData\Anaconda3\lib\site-packages\pySDC\implementations\sweeper_classes\boris_2nd_order.py", line 163, in integrate
f = P.build_f(L.f[j], L.u[j], L.time + L.dt * self.coll.nodes[j - 1])

File "C:\ProgramData\Anaconda3\lib\site-packages\pySDC\playgrounds\Boris\spiraling_particle_ProblemClass.py", line 102, in build_f
(f.elec.values + np.cross(part.vel.values[:, 0], f.magn.values[:, 0]))

ValueError: could not broadcast input array from shape (3,3) into shape (3) `

Suggest to add an argument indicating the current time node at which the solution lives to the restrict and prolong function of the SpaceTransfer class.

This would require modification of the lines in the BaseTransfer class where these two functions are called.

Line 108 for example would need to be changed to

tmp_u = [self.space_transfer.restrict(F.u[0], F.sweep.coll.nodes[0])]

correct?

Existing transfer classes would need to be modified to match the new signature, but the argument would be set to *unused.

Hello!

when setting this parameter https://github.com/Parallel-in-Time/pySDC/blob/master/pySDC/playgrounds/Boris/penningtrap_playground.py#L41 to 0 or 1, I get the same result. Is it OK?

I also set level_params['restol'] to 0.

Thanks!

For some reason, test_Q_transfer_minimal[Equidistant] always seems to lead to a Runtime Warning in _polyint, see https://github.com/Parallel-in-Time/pySDC/runs/7786219453?check_suite_focus=true#step:5:134. This seems to happen each time, but I have no idea why.. @tlunet?

The function u_exact is used to set the initial conditions at time 0, but then for many problems without analytic solutions, it just returns the initial conditions when called at a later time without raising an error.
I suggest adding assert 'exact solution is only available at t=0 for this problem!' to all problems that don't provide some kind of exact solution.

Alternatively, Scipy gives rather good reference solutions which can be computed very quickly, as is done here.

Derive a class imex_sweeper_scalar from imex_sweeper, which contains the functions to generate matrix representations of the sweeper.

Upon instantiation, this subclass should check that the attached problem is of the correct type.

Is there an easy way to run MLSDC/PFASST without applying the tau-correction? I am wondering what the easiest way is to set it to zero.

Looks as if everything is now in the master branch ... I deleted the branch from the PR, not sure why there is still a sharpclaw branch around. I suggest to delete it?

Adopt the class providing standard integrators to that it relies on the interface defined by the problem classes.

Collect all standard integrators at one place and implement a testing framework for them.

In the BaseTransfer class lines 108-110, why has F.u M+1 entries but F.sweep.coll.nodes only M? I assume that the first entry in F.u is the initial value, correct? This would mean we have

F.u[0] --> F.status.time (start time of time step)
F.u[1] --> F.sweep.coll.nodes[0] (first quadrature node)
F.u[2] --> F.sweep.coll.nodes[1] (second quadrature node)

etc. Is my interpretation correct?

For some numbers, the tests using the new node generator leads to errors when testing the equidistant distribution. Looks like it happens only for some random numbers and is thus hard to reproduce. One instance of this can be seen here: https://github.com/Parallel-in-Time/pySDC/runs/7713874352?check_suite_focus=true#step:5:339

@tlunet: If you have an idea, go ahead and fix. If not, then I'll fix the random seed and ignore it.

This leads to the side-effect that two calls of run_pfasst use the same stats object.. there seems to be no way to remove the content of the object before starting the second call, so the statistics may get mixed up between the calls..

At the moment, it is possible to instantiate an object of type imex_1st_order by providing only the parameters collocation_class and num_nodes.

However, when calling e.g. the integrate function, an error is thrown because no problem or level is attached. This should be caught during instantiation, but I am not sure what is the proper way to that.

Running nosetests throws the following warning, which seems to indicate that M1 and M2 are not necessarily integers. This is odd.

tests.test_collocation.TestCollocation.test_1 ... pySDC/pySDC/Collocation.py:85: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
L = np.zeros([M1, M2])
/usr/not-backed-up/pySDC/pySDC/Collocation.py:87: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
Lp = np.zeros(M2)
/usr/not-backed-up/pySDC/pySDC/Collocation.py:97: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
Lp = M2 * (L[:, M1 - 1] - y * L[:, M2 - 1]) / (1 - y ** 2)
/usr/not-backed-up/pySDC/pySDC/Collocation.py:99: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
y = y0 - L[:, M2 - 1] / Lp
/usr/not-backed-up/pySDC/pySDC/Collocation.py:96: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
L[:, k] = ((2 * k - 1) * y * L[:, k - 1] - (k - 1) * L[:, k - 2]) / k

The residual in the IMEX acoustic-advection example is not decreasing: It remains roughly at unity. The solution, however, appears to be correct.

Copy the corresponding functions from the Matlab msdc_oo code.

I'm pretty sure they way interpolation and restriction in time (i.e. between collocation nodes) is currently done is not a great idea, see here. @tlunet: Would you mind having a look at it? Can your interpolation routine make this better?

We could try to speed up the setup phase(s) by allocating/creating big matrices e.g. in the Boussinesq example just once, before the levels are set up. Then, pass these matrices to the ProgramClass and just add a link here.

Do we need to make numbers in formulas with real value output into reals, i.e. e.g.
here?

In preparing for importing pyclaw, I found an issue where the code runs properly if I change 1 to 1. there but crash with a 'division by zero' if not. The reason is that dx is zero and not set properly.

Unfortunately, I was not able to reproduce this in the python2 branch.

The time is an argument for both eval_fexpl and eva_fimpl but not for solve_system. For consistency, it should also be passed to solve_system (also, I need it for my time-dependent BC).

The new adaptivity feature in the controller(s) should be showcased in a dedicated project, including tests and documentation.

• Create new project infrastructure (folder, test, docs)
• Showcase adaptive SDC with and without HotRod (incl. tests and doc)
• Showcase adaptive multi-step SDC with and without HotRod (incl. tests and doc)
• Do all this in the parallel controller as well (incl. test and doc)

TBC.

Hello!

I have a question about implicit dependency in the SDC iteration.
Could you please explain how c_k value is defined in pySDC code?
I expected that formula (48) (Winkel, 2015) should be used.

In pySDC c_k value is defined in another way by https://github.com/Parallel-in-Time/pySDC/blob/master/pySDC/implementations/sweeper_classes/boris_2nd_order.py#L108.
After that it is transmitted to boris_solver at line 133, 136

I think there is a trick but I will be grateful for some comments.

The BaseTransfer class contains loads of functionality that assumes two meshes between which data is transfered, e.g. here

This might not be the case e.g. when using physics based coarsening or also for particle methods. In order to have a clearer description of the interface, does it may be make sense to move all this functionality into a derived class so that BaseTransfer becomes essentially an abstract class?

Running a Conda environment in Linux terminal, all dependencies installed without issue. Running 'spiraling_particle_playground.py' throws the following error for the progressbar:

(pySDC) [---User and PC name dedacted--- Boris]$ python spiraling_particle_playground.py
controller - WARNING: classic controller is about to become deprecated, use multigrid controller instead
Traceback (most recent call last):
File "spiraling_particle_playground.py", line 110, in
uinit, stats, a0 = main(dt=0.01, Tend=1)
File "spiraling_particle_playground.py", line 68, in main
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
File "/home/csunix/mn12kms/.conda/envs/pySDC/lib/python3.7/site-packages/pySDC/implementations/controller_classes/allinclusive_classic_nonMPI.py", line 96, in run
self.hooks.pre_run(step=S, level_number=0)
File "/home/csunix/mn12kms/.conda/envs/pySDC/lib/python3.7/site-packages/pySDC/playgrounds/Boris/spiraling_particle_HookClass.py", line 30, in pre_run
self.bar_run = progressbar.ProgressBar(max_value=progressbar.UnknownLength)
TypeError: init() got an unexpected keyword argument 'max_value'

Just updated my fork repository (https://github.com/tlunet/pySDC) to the latest version, but then it runs some test to check a "gitlab mirror", and it does not pass (see error log in https://github.com/tlunet/pySDC/runs/5378850003?check_suite_focus=true).
Seems to be some authentication issue with the gitlab mirror repository from Helmotz Codebase repo ...

Should we try to fix it or just ignore ?

Dependabot couldn't authenticate with https://pypi.python.org/simple/.

You can provide authentication details in your Dependabot dashboard by clicking into the account menu (in the top right) and selecting 'Config variables'.

View the update logs.

E.g. lines 34 & 36 in trajectory_HookClass.py fail with the error

TypeError: __init__() got an unexpected keyword argument 'max_value'

when running vanderpol_playground.py.

According to here, the keyword now has to be maxval instead of max_value.

However, fixing this causes the code to fail in line 51 with the error

AttributeError: next_update

Apparently, this is also due to a change in the progressbar module: see here.

At the moment, the upwind stencil is correct for negative advection velocity, i.e. from right to left. For the Boussinesq example, we want positive advection velocity, i.e. from left to right. Fix the stencil accordingly.

We receive new u[0], but u[1] is not changed (e.g. in the prolong-routines).. the same is true for f[1] and uold, fold in the prolongation.

Dependabot can't evaluate your Python dependency files.

As a result, Dependabot couldn't check whether any of your dependencies are out-of-date.

The error Dependabot encountered was:

InstallationError("Invalid requirement: 'fenics=2018.1.0'\n= is not a valid operator. Did you mean == ?")

You can mention @dependabot in the comments below to contact the Dependabot team.

Dependabot couldn't authenticate with https://pypi.python.org/simple/.

You can provide authentication details in your Dependabot dashboard by clicking into the account menu (in the top right) and selecting 'Config variables'.

View the update logs.

Where would be a good place to add a script that can be run to generate the nodes, weights, Q and S matrix for a specific type and number of nodes?

In the constructor of Sweeper.py, shouldn't there be asserts to make sure the later used entries of params exist? It seems that num_nodes and collocation_class have to exist in every set of parameters, independent of the underlying sweeper.

Therefore, it would make sense to check for their existence in the constructor of the sweeper class?

In runconvergence.py in acoustic_1d_imex:

• run simulation until Tend=1.0 with Nsteps = 64, i.e. dt = 0.015625. In this case, the output from the last step reads 2015-12-02 13:42:04,252 - Process 0 on time 0.984375 at stage IT_FINE_SWEEP: Level: L0 -- Iteration: 1 -- Residual: 1.87173732e-03 2015-12-02 13:42:04,260 - Process 0 on time 0.984375 at stage IT_FINE_SWEEP: Level: L0 -- Iteration: 2 -- Residual: 3.94070291e-05 2015-12-02 13:42:04,270 - Process 0 on time 0.984375 at stage IT_FINE_SWEEP: Level: L0 -- Iteration: 3 -- Residual: 8.01802483e-07 2015-12-02 13:42:04,279 - Process 0 on time 0.984375 at stage IT_FINE_SWEEP: Level: L0 -- Iteration: 4 -- Residual: 1.57838695e-08. The time matches Tend-dt = 1.0 - 0.015625 = 0.984375
• when chaning to Nsteps=65, the output from the last step becomes: 2015-12-02 13:47:50,509 - Process 0 on time 1.000000 at stage IT_FINE_SWEEP: Level: L0 -- Iteration: 1 -- Residual: 1.80955876e-03 2015-12-02 13:47:50,518 - Process 0 on time 1.000000 at stage IT_FINE_SWEEP: Level: L0 -- Iteration: 2 -- Residual: 3.72419523e-05 2015-12-02 13:47:50,527 - Process 0 on time 1.000000 at stage IT_FINE_SWEEP: Level: L0 -- Iteration: 3 -- Residual: 7.40115315e-07 2015-12-02 13:47:50,536 - Process 0 on time 1.000000 at stage IT_FINE_SWEEP: Level: L0 -- Iteration: 4 -- Residual: 1.42181506e-08 and the plotted solution suddenly deviates significantly (i.e. visibly in the generated plot) from the exact solution. For some reason, the final time is now exactly Tend instead of Tend-dt as in the case above!

The only difference in parameters in both simulations is the increase in time steps from 64 to 65.

To ease the setup and continuous usage of pySDC, I recommend using a requirements.txt file, which can be directly used by pip to install all dependencies.