Not all particle properties are initialized correctly

[bangerth]

We initialize particle properties based on fields at the very beginning of the simulation via

signals.post_set_initial_state
World::setup_initial_state ()
World::initialize_particles()
World::local_initialize_particles()
PropertyManager::initialize_one_particle()

and the signal is triggered here:

  start_time_iteration:

    if (parameters.resume_computation == false)
      {
        TimerOutput::Scope timer (computing_timer, "Setup initial conditions");

        timestep_number           = 0;
        time_step = old_time_step = 0;

        if (! parameters.skip_setup_initial_conditions_on_initial_refinement
            ||
            ! (pre_refinement_step < parameters.initial_adaptive_refinement))
          {
            // Add topography to box models after all initial refinement
            // is completed.
            if (pre_refinement_step == parameters.initial_adaptive_refinement)
              signals.pre_set_initial_state (triangulation);

            set_initial_temperature_and_compositional_fields ();
            compute_initial_pressure_field ();

            signals.post_set_initial_state (*this);
          }
      }

    // Start the principal loop over time steps. At this point, everything
    // is completely initialized, so set that status as well
    simulator_is_past_initialization = true;
    do
     {
        [main time loop]

The issue is that at that time, not all fields are already initialized. Specifically, we don't have a pressure field that makes sense unless the compute_initial_pressure_field() function computes anything sensible.

What does this matter? If one initializes particle properties from these fields, these fields can be the initial conditions for ODEs and one would imagine that we don't quite get the right solution of that ODE if we have the wrong initial conditions. An example is the "pT" particle property that, for problems without a hydrostatic pressure, are simply initialized with a zero pressure.

What to do: Would it make sense to initialize particle properties after the first Stokes solution step? I think that's wrong too since (i) we solve for compositional fields before the Stokes system, and (ii) I believe that particle properties can enter the compositional fields. Ideas?

(This all was pointed out to me by @cedrict .)