Can libmesh solve several equation systems in parallel?

[nuomi68]

[nuomi68]

Could you print out solutions u and v respectively in examples/systems_of_equations/ex4?

[roystgnr]

The code you ... screenshotted, for some reason? ... solves in serial, not parallel. To solve coupled systems of equations in parallel, they need to be in the same System. To solve independent systems of equations in parallel, you can put them in separate EquationSystems, each with a different subcommunicator of the same global communicator, and solve each on the processors assigned to that subcommunicator. (this is sometimes useful for Monte Carlo or other "embarrassingly parallel" statistics methods).

[roystgnr]

I'm afraid I don't. I see a couple unit tests of Communicator::split() that I think were added around that time, but not even unit tests of solves with the resulting subcommunicators, and it looks like the apps that were doing that didn't get open sourced.

Which isn't a great sign, since they're not in CI either. In theory, you simply split() a communicator, you initialize your Mesh with the subcommunicator instead of the full communicator, and then the rest of your code looks like ordinary libMesh code, except that where you have any custom communication code you have to be careful about where you use the subcommunicator (for the parallel solves) vs the full communicator (to aggregate your results afterward). But in practice, if that's not covered in CI there might have been some regression we missed somewhere along the line (though I can't imagine what off the top of my head).

[jwpeterson]

At the risk of getting too far ahead of ourselves, can you give us a bit more information about what you are actually trying to do? Solving different Systems at the same time in parallel isn't a standard libMesh operation, and it isn't going to be trivial to do (@roystgnr's nice explanation notwithstanding).

[nuomi68]

I'm happy to share our thoughts, which may be unrealistic. We intend to use libmesh to solve the neutron transport equation. For the SN discrete method, the equation needs to be solved once for each Angle, and it may be possible to solve the equation for each Angle in parallel when the coupling is ignored. Maybe our idea is not good, or maybe we don't need multiple systems in parallel, and we can solve it in parallel within the same system.

[roystgnr]

That's a clever idea. You definitely can solve it in parallel within the same system (perhaps using a custom CouplingMatrix to make sure the zero-coupling blocks are gone, not just full of 0.0), and I'm actually not sure what would be more efficient. If you set up a separate System for each angle, then you risk one kind of load imbalance (some subcommunicators stuck twiddling their thumbs because their solve converged faster than the others') that's hard to fix; if you set up a single System then you risk a different kind (each iteration giving equal time to every DoF even though some don't need it) that's impossible to fix. The single System code can be written require fewer flops (you're doing much of the same calculations for every angle and don't need to repeat them), and it'll probably take less RAM (no need for multiple copies of your mesh), but it'll involve more communication (you'll need a partition into fewer elements per processor, assuming the same total number of processors).

I'd write a single System code to start with, just because that'll be easier. If you do try the embarrassingly-parallel angle, though, I'd love to know what the results look like.

[nuomi68]

Thank you so much for your insightful response! I'm inclined to follow your advice and start with a single System code due to its potential for reduced flops and memory usage. And may be easy to implement, relative to multiple systems. Also, using a single system rather than multiple systems is definitely better for initial implementation. If successful, I'll be more than happy to share the results with you .

[nuomi68]

Our code is as follows, currently considering a very simplified neutron transport equation, so when solving each Angle, only need to change the discrete Angle direction and weight, the other values do not change, so the execution of num_N loop parallel is possible.

  std::vector<std::vector<Real>> Omiga = {{0.350021,0.350021,0.868890}, {0.350021,0.868890,0.350021},
                                            {0.868890,0.350021,0.320021}};
  std::vector<Real> Weight = {1.0/3.0,1.0/3.0,1.0/3.0};
for(auto iter : make_range(origin_iters))
  {                 
 
    phi_bar.operator=(*bar_cache.clone());
    bar_cache.zero();

    for(auto N : make_range(num_N))
    {

      equation_systems.parameters.set<std::vector<Real>>("omiga")  = Omiga[N];
      equation_systems.parameters.set<Real>("weight") = Weight[N];

      equation_systems.get_system("neutron transport").solve();

      NumericVector<Number> &  solution=(*system.current_local_solution); 
      \\calculate the mean neutron flux, prepare for the next outer iteration
      bar_cache.operator+=(solution.operator*=(Weight[N]));
    }
  }
}