Allow external initialization of MPI and p4est

[sloede]

This PR allows some packages, namely MPI and p4est, to be initialized (and finalized) externally. That is, in the constructor for Utilities::MPI_InitFinalize, either library is only initialized if it has not already been initialized before. For consistency, if external initialization was detected, finalization in the destructor is also omitted.

The motivation for this proposed change is a joint project with @kronbichler and @ranocha, where we use deal.II as a library and not as the main driver of a simulation. Since MPI and p4est are already initialized when deal.II is loaded, we need to be able to switch off the automatic initialization of third-party libraries.

An optional function parameter is added to the constructor of Utilities::MPI_InitFinalize that enables the new behavior. By default, everything remains as is, i.e., for existing codes nothing should change and, e.g., an already initialized MPI library will cause an error to be thrown as before.

Please note that I am not an expert on the inner workings of deal.II, and that thus there might be smarter ways to achieve the same effect. Therefore, please view this PR not necessarily as the best possible implementation but rather as a means to get the discussion started.

[kronbichler]

One thing I would like to add to the discussion is the fact that @sloede and @ranocha have first thought we might completely skip using MPI_InitFinalize() but then discovered that deal.II has several additional facilities located in MPI_InitFinalize(), e.g.

dealii/source/base/mpi.cc

Lines 780 to 828 in 03e8a82

	// Initialize Kokkos
	{
	// argv has argc+1 elements and the last one is a nullptr. For appending
	// one element we thus create a new argv by copying the first argc
	// elements, append the new option, and then a nullptr.
	std::vector<char *> argv_new(argc + 2);
	for (int i = 0; i < argc; ++i)
	argv_new[i] = argv[i];
	std::stringstream threads_flag;
	#if KOKKOS_VERSION >= 30700
	threads_flag << "--kokkos-num-threads=" << MultithreadInfo::n_threads();
	#else
	threads_flag << "--kokkos-threads=" << MultithreadInfo::n_threads();
	#endif
	std::string threads_flag_string = threads_flag.str();
	argv_new[argc] = const_cast<char *>(threads_flag_string.c_str());
	argv_new[argc + 1] = nullptr;
	// The first argument in Kokkos::initialzie is of type int&. Hence, we
	// need to define a new variable to pass to it (insted of using argc+1
	// inline).
	int argc_new = argc + 1;
	Kokkos::initialize(argc_new, argv_new.data());
	}
	// we are allowed to call MPI_Init ourselves and PETScInitialize will
	// detect this. This allows us to use MPI_Init_thread instead.
	#ifdef DEAL_II_WITH_PETSC
	# ifdef DEAL_II_WITH_SLEPC
	// Initialize SLEPc (with PETSc):
	finalize_petscslepc = SlepcInitializeCalled ? false : true;
	ierr = SlepcInitialize(&argc, &argv, nullptr, nullptr);
	AssertThrow(ierr == 0, SLEPcWrappers::SolverBase::ExcSLEPcError(ierr));
	# else
	// or just initialize PETSc alone:
	finalize_petscslepc = PetscInitializeCalled ? false : true;
	ierr = PetscInitialize(&argc, &argv, nullptr, nullptr);
	AssertThrow(ierr == 0, ExcPETScError(ierr));
	# endif
	// Disable PETSc exception handling. This just prints a large wall
	// of text that is not particularly helpful for what we do:
	PetscPopSignalHandler();
	#endif
	// Initialize zoltan
	#ifdef DEAL_II_TRILINOS_WITH_ZOLTAN
	float version;
	Zoltan_Initialize(argc, argv, &version);
	#endif

or

dealii/source/base/mpi.cc

Lines 963 to 999 in 03e8a82

	// Before exiting, wait for nonblocking communication to complete:
	for (auto request : requests)
	{
	const int ierr = MPI_Wait(request, MPI_STATUS_IGNORE);
	AssertThrowMPI(ierr);
	}
	// Start with deal.II MPI vectors and delete vectors from the pools:
	GrowingVectorMemory<
	LinearAlgebra::distributed::Vector<double>>::release_unused_memory();
	GrowingVectorMemory<LinearAlgebra::distributed::BlockVector<double>>::
	release_unused_memory();
	GrowingVectorMemory<
	LinearAlgebra::distributed::Vector<float>>::release_unused_memory();
	GrowingVectorMemory<LinearAlgebra::distributed::BlockVector<float>>::
	release_unused_memory();
	// Next with Trilinos:
	# ifdef DEAL_II_WITH_TRILINOS
	GrowingVectorMemory<
	TrilinosWrappers::MPI::Vector>::release_unused_memory();
	GrowingVectorMemory<
	TrilinosWrappers::MPI::BlockVector>::release_unused_memory();
	# endif
	#endif
	// Now deal with PETSc (with or without MPI). Only delete the vectors if
	// finalize hasn't been called yet, otherwise this will lead to errors.
	#ifdef DEAL_II_WITH_PETSC
	if (!PetscFinalizeCalled)
	{
	GrowingVectorMemory<
	PETScWrappers::MPI::Vector>::release_unused_memory();
	GrowingVectorMemory<
	PETScWrappers::MPI::BlockVector>::release_unused_memory();
	}

. So I guess some infrastructure is needed to cope with this use case (which I think is important as a library). But similar things must have been done in the xSDK setup - how is the initialization handled there?

[Rombur]

The goal of xSDK is to be able to compile all the libraries together. That's it.

Another solution would be to create a new class that only initializes/finalizes deal.II and let the user take care of initializing all the libraries they need. The advantage is that the code would be more clean (we don't have if everywhere). With this PR, we will have a class called MPI_InitFinalize that does everything but initializing MPI which I find awkward. I am not opposed to this PR but it's not my favorite solution.

[masterleinad]

But similar things must have been done in the xSDK setup - how is the initialization handled there?

xSDK just makes sure that deal.II can be built with other xSDK packages like PETSc and Trilinos. There is no application/test case that uses deal.II as a library.

[kronbichler]

With this PR, we will have a class called MPI_InitFinalize that does everything but initializing MPI which I find awkward. I am not opposed to this PR but it's not my favorite solution.

I agree, that would be awkward. So that makes me lean towards a different class/struct to collect the infrastructure needed by deal.II proper (maybe including Kokkos, maybe not), and let everything else like PETSc, p4est, MPI be handled by the user?

[tjhei]

So that makes me lean towards a different class/struct to collect the infrastructure needed by deal.II proper

Agreed. Something like dealii::InitFinalize() that can be added as a member to MPI_InitFinalize.

[Rombur]

@kronbichler yes, I like it better.

I would let the user initialize Kokkos. When using the bundled version of Kokkos, we only enable the serial backend. So not only we don't need the code that Daniel added but you don't need to initialize Kokkos yourself (we safely do it for you). It's only when using OpenMP or a GPU backend that you need to do something special. If users are compiling their own Kokkos, we can let them initialize Kokkos themselves.

[masterleinad]

I agree that the new initialization function shouldn't initialize Kokkos.

[sloede]

I've followed the discussion here with great interest.

So that makes me lean towards a different class/struct to collect the infrastructure needed by deal.II proper (maybe including Kokkos, maybe not), and let everything else like PETSc, p4est, MPI be handled by the user?

From a research software engineering perspective, this makes sense to me too. However, I have to admit that I am completely out of my depth with the proposed new approach, thus I am afraid I will not be able to implement this myself. I'd be happy to test drive an implementation from someone else for our application though, to give a user's perspective.

[bangerth]

How about a class dealii::InitFinalize that takes a list of things in the constructor (an enum's elements, merged via operator|) that need to be initialized? Or maybe better for the current context, with a list of base libraries that should not be initialized? This class basically does what the current one does, except that one can in/exclude individual code blocks.

Then the old MPI_InitFinalize name could be a class derived from that new class, calling the base class's constructor with "all dependency libraries should be initialized".

[stefanozampini]

I don't understand why you need to require an extra flag allow_external_initialization. For example, the PETSc/SLEPC part can be initalized outside of this function and deal.ii does not complain. What is different for MPI?

[sloede]

Currently in master, deal.II does not allow MPI to be initialized outside of deal.II:

dealii/source/base/mpi.cc

Lines 753 to 757 in 3c37cfe

	int MPI_has_been_started = 0;
	ierr = MPI_Initialized(&MPI_has_been_started);
	AssertThrowMPI(ierr);
	AssertThrow(MPI_has_been_started == 0,
	ExcMessage("MPI error. You can only start MPI once!"));

I thus assumed that the design goals for this function were that any abnormal usage by the user should result in a fast and hard failure, even though in this case it would have been possible to, e.g., first check the current thread level and see if it is sufficient for deal.II. This makes sense to me - I have used the same approach many times before, assuming that if anything is happening out of the ordinary, the user should probably be informed immediately instead of potentially running an expensive simulation and failing after a couple of hours.

My approach in this PR was therefore motivated by wanting to keep this strict safety against user errors as the default. I therefore made them set a flag explicitly, which indicates to deal.II "I know what I'm doing" but also provides the feedback to the users "if you had to set a flag for this, you're probably in dragon land now".

However, I am not a regular committer to deal.II, so everything I wrote was based only on assumptions. I'd be more than happy to make this work without a flag, if this is so desired. However, from the convo in this PR, I already got the feeling that a very different approach is to be followed for allowing external initialization, so I don't know if this question will not be rendered moot anyway.

[stefanozampini]

same here

[bangerth]

How do we want to proceed with this PR? Did #15572 address at least some parts of what you wanted to do?

[sloede]

How do we want to proceed with this PR? Did #15572 address at least some parts of what you wanted to do?

Thanks for your effort! I have not yet found the time to review and understand how I can do with #15572 what we originally tried to achieve, but it is on my TODO list 😊