Coupling with load balancing causes noticeable MaMiCo overhead in sparse scenarios

[amartyads]

When using static load balancing with ls1+autopas, the force calculation time reduces, but the time taken by the mamico steps per MD timestep increases. It seems like mamico only cares about number of coupling cells per rank and not number of molecules per rank, and thus, in a state balanced by molecules, the cells are not balanced. So if the domain size is small, then the mamico runtime dominates, and loadbalanced scenarios are performing worse.

This was noticed while trying to balance a droplet scenario. For a couette scenario, this wasn't noticed. Therefore, for denser setups, this isn't an issue, since the coupling overhead is sufficiently smaller than the computation.

Some example outputs (all mamico+ls1+autopas, 2 md instances, and my multiphase droplet scenario, and mamico cellsize 5) are given below. Compare the force calculation time to the MamicoCoupling plugin time. The computation time combines both, despite the weird alignment. All particle counts are estimated. All temperatures are 0.8.

For MD120, relevant values:
Droplet diameter: 80
Liquid density: 0.73056664
Vapour density: 0.01991599
MD120, 4 ranks per axis, equally partitioned, 128 ranks,
per subdomain: max particles: 18213 min particles: 538 max density: 0.6745469705 min density: 0.01991599000

INFO:   2025-03-30T00:57:32 161.93   [0]        Computation in main loop took: 80.3111 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                Decomposition took: 5.47365 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                        Container update took: 0.254231 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                        Communication took: 5.21271 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                                initSend() took: 0.442844 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                                testRecv() took: 0.0856563 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                        Cache update took: 0.000861809 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                Computation took: 74.807 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                        Force calculation took: 26.7712 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                IO in main loop took: 0.0239461 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                        Plugin MamicoCoupling took: 24.1269 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                Enforcing boundary conditions took: 0.00191764 sec
INFO:   2025-03-30T00:57:32 161.93   [0]        Timer SIMULATION_IO took: 0 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                Phasespace creation took: 0.306877 sec
INFO:   2025-03-30T00:57:32 161.93   [0]                Final IO took: 1.607e-05 sec

MD120, 4 ranks per axis, balanced in ratio x="9 ; 3 ; 3 ; 9" y="9 ; 3 ; 3 ; 9" z="9 ; 3 ; 3 ; 9", 128 ranks
per subdomain: max particles: 3889 min particles: 2466 max density: 0.7305666400 min density: 0.03491958966

INFO:   2025-03-30T00:50:55 296.865  [0]        Computation in main loop took: 149.326 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                Decomposition took: 2.80304 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                        Container update took: 0.793163 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                        Communication took: 2.00244 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                                initSend() took: 1.31857 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                                testRecv() took: 0.489584 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                        Cache update took: 0.000848305 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                Computation took: 146.464 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                        Force calculation took: 6.32611 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                IO in main loop took: 0.0528344 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                        Plugin MamicoCoupling took: 139.689 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                Enforcing boundary conditions took: 0.0026409 sec
INFO:   2025-03-30T00:50:55 296.865  [0]        Timer SIMULATION_IO took: 0 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                Phasespace creation took: 0.29144 sec
INFO:   2025-03-30T00:50:55 296.865  [0]                Final IO took: 1.7023e-05 sec

For MD180, relevant values:
Droplet diameter: 100
Liquid density: 0.73049450
Vapour density: 0.01988934
md180, 4 ranks per axis, equally partitioned, 64 ranks
per subdomain: max particles: 46300 min particles: 1812 max density: 0.5080906812 min density: 0.01988934000

INFO:   2025-03-30T01:52:45 264.836  [0]        Computation in main loop took: 253.426 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                Decomposition took: 5.81495 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                        Container update took: 0.657443 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                        Communication took: 5.15084 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                                initSend() took: 1.14164 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                                testRecv() took: 0.148227 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                        Cache update took: 0.000795103 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                Computation took: 247.58 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                        Force calculation took: 60.3213 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                IO in main loop took: 0.025381 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                        Plugin MamicoCoupling took: 126.98 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                Enforcing boundary conditions took: 0.00183197 sec
INFO:   2025-03-30T01:52:45 264.836  [0]        Timer SIMULATION_IO took: 0 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                Phasespace creation took: 0.641713 sec
INFO:   2025-03-30T01:52:45 264.836  [0]                Final IO took: 1.9647e-05 sec

md180, 4 raNks per axis, ratio 10,8,8,10 for x, 9,9,9,9 for y and z, 64 ranks
per subdomain: max particles: 44169 min particles: 1611 max density: 0.5452963200 min density: 0.01988934000

INFO:   2025-03-30T01:58:04 254.864  [0]        Computation in main loop took: 243.041 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                Decomposition took: 2.29667 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                        Container update took: 0.735421 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                        Communication took: 1.55442 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                                initSend() took: 1.27205 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                                testRecv() took: 0.158925 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                        Cache update took: 0.000798889 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                Computation took: 240.711 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                        Force calculation took: 57.7005 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                IO in main loop took: 0.0268437 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                        Plugin MamicoCoupling took: 162.466 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                Enforcing boundary conditions took: 0.00200887 sec
INFO:   2025-03-30T01:58:04 254.864  [0]        Timer SIMULATION_IO took: 0 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                Phasespace creation took: 0.688376 sec
INFO:   2025-03-30T01:58:04 254.864  [0]                Final IO took: 2.1109e-05 sec

md180, 4 ranks per axis, ratio 10,8,8,10 for all x, y and z, 64 ranks
per subdomain: max particles: 39972 min particles: 1989 max density: 0.6245643667 min density: 0.01988934000

INFO:   2025-03-30T03:36:57 293.6    [0]        Computation in main loop took: 281.852 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                Decomposition took: 2.77866 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                        Container update took: 0.901442 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                        Communication took: 1.87055 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                                initSend() took: 1.56204 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                                testRecv() took: 0.183096 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                        Cache update took: 0.000794358 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                Computation took: 279.04 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                        Force calculation took: 54.1327 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                IO in main loop took: 0.0267248 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                        Plugin MamicoCoupling took: 224.43 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                Enforcing boundary conditions took: 0.0018465 sec
INFO:   2025-03-30T03:36:57 293.6    [0]        Timer SIMULATION_IO took: 0 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                Phasespace creation took: 0.689739 sec
INFO:   2025-03-30T03:36:57 293.6    [0]                Final IO took: 1.8715e-05 sec

[Thinkpiet]

Can you post some config files, to give more info and make this reproducible?

[amartyads]

I think the easiest way to replicate would be to

1. Checkout my repo https://github.com/amartyads/singleDropletBuilder
2. Edit config.json appropriately. Accepted jobmanagers are pbs slurm, accepted systems are hawk hsuper
3. Run create-droplet-config.py
4. Change <loglevel> in coupled/ls1config.xml to INFO
5. Run add-and-run-jobscripts.py. I haven't comprehensively tested the jobscripts and stuff on hsuper, but it should be fine.

[amartyads]

More outputs after Piet's suggestions:
Same md180 as above (final config with balancing 10,8,8,10 all ranks), with the CFD split into 2;2;2 ranks, and MD droplet centered in the CFD domain (it already was):

INFO:   2025-04-01T18:41:23 278.13   [0]        Computation in main loop took: 276.792 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                Decomposition took: 2.73312 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                        Container update took: 0.886147 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                        Communication took: 1.84018 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                                initSend() took: 1.53164 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                                testRecv() took: 0.176891 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                        Cache update took: 0.000809443 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                Computation took: 274.025 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                        Force calculation took: 53.0979 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                IO in main loop took: 0.0277559 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                        Plugin MamicoCoupling took: 220.448 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                Enforcing boundary conditions took: 0.00215461 sec
INFO:   2025-04-01T18:41:23 278.13   [0]        Timer SIMULATION_IO took: 0 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                Phasespace creation took: 0.645489 sec
INFO:   2025-04-01T18:41:23 278.13   [0]                Final IO took: 1.7964e-05 sec

Same md180 as above (final config with balancing 10,8,8,10 all ranks), with only 1 MD instance and 64 ranks:

INFO:   2025-04-01T18:52:10 282.273  [0]        Computation in main loop took: 276.042 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                Decomposition took: 2.70963 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                        Container update took: 0.879042 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                        Communication took: 1.82396 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                                initSend() took: 1.52373 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                                testRecv() took: 0.176409 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                        Cache update took: 0.000787417 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                Computation took: 273.3 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                        Force calculation took: 53.2086 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                IO in main loop took: 0.0260154 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                        Plugin MamicoCoupling took: 219.629 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                Enforcing boundary conditions took: 0.00196274 sec
INFO:   2025-04-01T18:52:10 282.273  [0]        Timer SIMULATION_IO took: 0 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                Phasespace creation took: 0.62351 sec
INFO:   2025-04-01T18:52:10 282.273  [0]                Final IO took: 1.7773e-05 sec

[amartyads]

Summary of chat with Piet regd this issue:

1. Added densities and particle counts, added results of running parallel CFD and single MD instance (no changes seen)
2. Mamico output is turned off, no output other than stdout
3. It would be nice to see the mamico coupling overhead per MD timestep, to see if it's consistent or if there are spikes
4. Mamico does some operations per cell. And because mamico needs to traverse all the cells. That has never been very efficient, they can be scattered in memory and are traversed by pointers. For a reasonable number of particles per cell it does not matter, since then time for traversing particles and force calculation is more heavy. But for small number of particles per cell one would expect to see the effects of cell traversal and cell wise operations.
5. We haven't tested with sparse gases before
6. We need benchmarks for cell traversals and cell mapping on mamico side
7. We should write a benchmark that creates a very large grid of coupling cells in macroscopic cell service, then puts exactly zero particles in each grid cell, then calls every cell mapping (ca 15 different ones) for every coupling cell, then outputs the run time. And then we should try to improve mamico so that that time goes towards zero.

I will find the overhead per MD timestep later, and update the issue.

[amartyads]

I timed the MamicoCoupling plugin to check if there's one large spike or whether the slowdown is cumulative, and it is indeed cumulative. Full output here:

output-108-3D-1MD-lotout.txt

And here is output from all ranks (much larger file):

output-108-3D-1MD-lotout.txt

The config is same md180 as above (final config with balancing 10,8,8,10 all ranks), with only 1 MD instance and 64 ranks.

Mamico coupling happens in two phases, i.e. immediately before and after force calculations. This is what happens:

void MamicoCoupling::beforeForces(ParticleContainer *particleContainer, DomainDecompBase *domainDecomp,
								  unsigned long simstep) {
	if (_couplingEnabled) {
		// This object should be set by MaMiCo after the plugins are created in the simulation readxml file.
		// Even though this method is called before the object is set, at this point the coupling switch is always off
		_couplingCellService->processInnerCouplingCellAfterMDTimestep();
		_couplingCellService->distributeMass(simstep);
		_couplingCellService->applyTemperatureToMolecules(simstep);
#ifndef MARDYN_AUTOPAS
		particleContainer->deleteOuterParticles();
#endif
	}
}

void MamicoCoupling::afterForces(ParticleContainer *particleContainer, DomainDecompBase *domainDecomp,
								 unsigned long simstep) {
	if (_couplingEnabled) {
		_couplingCellService->distributeMomentum(simstep);
		_couplingCellService->applyBoundaryForce(simstep);
	}
}

I've timed the beforeforces and afterforces separately every timestep, and at the end I have the cumulatives:

INFO:   2025-04-02T12:06:47 291.045  [0]        Cumulative before: 76.3575
INFO:   2025-04-02T12:06:47 291.045  [0]        Cumulative after: 149.401
INFO:   2025-04-02T12:06:47 291.045  [0]        Max before: 0.0886377
INFO:   2025-04-02T12:06:47 291.045  [0]        Max after: 0.172808

There are no spikes, it's a steady slowdown.

The [0] in the output line denotes rank. Some ranks are consistently slower. To translate rank to 3d coordinate, topology should not matter, since everything is symmetrical. Rank 0 for instance has 10x10x10 cells, and consistently takes longer than some other ranks, despite having fewer particles.

[amartyads]

Update: I tried a scenario with droplet size 80, with and without AutoPas. The problem seems to exist only if AutoPas is enabled.

Configurations are in 3D, for every axis. They denote ratios, just like above.
Specifics:
Done on HSUper, 64 ranks, 4x4x4 breakdown, periodic, no insertions/boundary forces, thermostat only on outer 2 layers
drop size 80, MD size 120, wall height 140

With autopas:

Configuration	Walltime(s)	ForceComputation(s)	CouplingPlugin(s)
default	277.511	66.8116	24.6058
7-5-5-7	190.224	46.4425	45.6804
8-4-4-8	219.087	25.4479	81.1048
9-3-3-9	332.866	16.3469	144.4670
10-2-2-10	557.091	21.0193	246.6410

Without autopas:

Configuration	Walltime(s)	ForceComputation(s)	CouplingPlugin(s)
default	186.3170	57.7266	4.12144
7-5-5-7	127.4600	39.4247	6.87380
8-4-4-8	79.4758	19.8167	10.17680
9-3-3-9	150.2280	58.3286	6.52194
10-2-2-10	144.7220	57.7209	7.74836

Thus more investigation needs to done on the coupling side. If I find evidence that the issue is due to a bug/lack of feature on AutoPas side, I'll close this issue.