Hi Wolfgang,
Apologies for the late reply. Rene had mentioned that is not immediately possible but not impossible in principle. I agree it might be clearer to you if I explain the problem I am trying to solve and why I asked about the Stokes’ solvers.
Main idea: I am basically looking at how the dynamic topography at the coremantle boundary (CMB) influences global mantle convection through the gravitational collapse of a rheologically weak layer.
The details: Buoyancydriven flows can create dynamic topography at the CMB where mantle downwellings “push" material slightly deeper into the core compared to the ambient mantle to create topographic lows. We believe that at these lows, mantle material will physically interact with the liquid metal to create a rheologically weak layer (with a density and dynamic viscosity that is in between the lower mantle and outer core) that is susceptible to lateral gravitational collapse. However, due to the persistence of large scale mantle buoyancy forces that vary on much longer time scales, there is a tendency to restore dynamic topography. Hence, lateral spreading of this weak layer would be able to draw additional rock downwards from the mantle. We also believe that the lateral spreading creates a secondary flow field that is superimposed on the buoyancy flow field to advect whatever scalar field(s) in the ambient mantle.
To solve this problem, we thought of the following steps during each iteration:

Obtain the first flow field v due to density variations in the mantle under a certain boundary condition e.g. free slip.

Calculate the dynamic topography along the CMB e.g. assuming the lower boundary is stressfree.

From the dynamic topography, obtain some velocity values along the CMB.

Use these velocity values as BCs to solve the Stokes’ equation to get the second flow field u .

Advect the relevant scalar fields (such as temperature) using u + v .

Repeat.
So in each iteration, the Stokes’ solver is called twice in steps 1 and 4 with different BCs. I have already done this in Matlab for a Boussinesq fluid that is isoviscous in 2D by solving the associated Poisson and Laplace equations using the vorticitystreamfunction formulation. The reason my supervisor and I want to use Aspect is because of the grid refinement capabilities that are available. We wish to resolve details at the boundary layer close to the CMB without needing to do the same for the ambient mantle as much.
I hope this gives you a clearer background to the question.
Best,
Kang Wei