I would like to find out if the following workflow I have in mind is feasible to implement in Aspect at the moment. So for every time step, this basically happens:
Solve the Stokes’ equations with one set of boundary conditions (BCs) over the whole domain to obtain one set of velocity solutions (let’s call this A).
Use A to perform some calculations to get a second set of BCs.
Solve the Stokes’ equations again but with the abovementioned BCs to obtain a second set of velocity solutions (let’s call this B).
Combine A and B to advect the relevant fields to the next time step.
I am interested to find out if such a scheme is supported at the moment by a plugin for instance or will it likely to involve heavy tinkering to the simulator if it is not.
Kang Wei – I think that’s a difficult setup to make work in ASPECT without a lot of internal modifications. You’d in essence have to fork the ASPECT code base to make that work.
But, my experience is that it’s worth starting not with how you want to do something, but with what it is and why. Can you explain what problem you are trying to solve and why it requires two Stokes solves? Maybe we can find a different way to achieve what you are looking for, before we think about how to implement it!
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 core-mantle boundary (CMB) influences global mantle convection through the gravitational collapse of a rheologically weak layer.
The details: Buoyancy-driven 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 stress-free.
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 vorticity-streamfunction 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.
I would think this could be implemented as a new NonlinearSolver scheme, but it requires several modifications inside ASPECT.
I am not sure how familiar you are with the internals of ASPECT and how boundary conditions work, though. This certainly requires understanding it and doing some implementation work.
Hi Kang,
I’m sure one could shoehorn this into ASPECT in the way Timo suggests. But I’m not convinced that that is the right approach, primarily because I don’t quite understand yet why you need the two solves to begin with.
Can you explain a bit more how your steps 3&4 look like? How do you determine the velocity field along the CMB boundary, and why could this velocity not be determined before the first solve and then used as the velocity boundary condition for step 1 instead? In other words, instead of decomposing the velocity field into two components, each of which is driven by a different physical effect, why not compute just one velocity field that has both effects as right hand side/boundary condition?
Maybe separately, I often find it useful to actually write out the equations one wants to solve. It often clarifies what needs to actually happen. In your case, could you write out what the equations are that describe the CMB velocity field?
I suppose that first solve is required to obtain the dynamic topography h given an initial buoyancy field of the mantle. Once the dynamic topography is obtained, the rate of change of h (due the the gravitational collapse of the rheologically weak material residing in the topographic lows h<0) gives the boundary condition of u_z at z=0. I will rewrite steps 1 to 5 with equations assuming an incompressible and isoviscous mantle that behaves as a Boussinesq fluid:
Solve \eta\nabla^2 \textbf{v}-\nabla P_\textbf{v} = \rho_0[1-\alpha(T-T_0)]g\hat{\textbf{z}} and \nabla\cdot\textbf{v} = 0 with free-slip boundary conditions
Calculate the dynamic topography h = \frac{\sigma_{zz}(z=0)}{\Delta\rho g} (h>0 denotes upward deformation)
Let u_z(x,y,z=0) = u_{z,0}(x,y). To obtain u_{z,0},
K is a constant and \nabla^2_H = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}.
Solve \eta\nabla^2\textbf{u} - \nabla P_\textbf{u} = 0 and \nabla \cdot \textbf{u} = 0 with free-slip boundaries at the top and sides. For the bottom boundary, u_z will be the value calculated in the previous step, while u_x and u_y will be shear-stress free.
To advect T, the transport equation is \frac{\partial T}{\partial t} + (\textbf{u}+\textbf{v})\cdot\nabla T= \kappa\nabla^2T
The steps above outline one iteration. I guess it is possible to just solve the Stokes equation once by moving steps 2 and 3 above step 1, and then just solve for the overall velocity field, but that would mean the dynamic topography values have to be obtained prior to the solve. Is that possible?
Kang,
yes, that is in fact how we typically do things: You solve only once and the velocity so computed is the result of (i) the buoyancy forces based on the temperature field computed in a previous step and (ii) the current shape of the domain, which has also been computed in a previous step. Previous step could be a temperature solve at the beginning of the time step computed before we get around to solving for the velocity field, or the domain shape update that happens at the end of the previous time step. In any case, the velocity so computed takes into account all effects that might cause motion, not just buoyancy.
I think you should take a look at the free surface module in ASPECT. It is used in several of the cookbooks, and extensively documented in a paper by Timo Heister and Ian Rose.
I’d like to point out one of the other reasons why I think that solving only once is the right way to go. Right now, you’re only using a constant viscosity that does not depend on the strain rate. But the real earth has a strain-rate dependent rheology. The problem is that you can only compute the viscosity if you know the velocity field, but in your case neither u not v alone corresponds to the velocity – it’s only u+v, which is not actually available in either of the two sub-steps you propose. You could make an assumption that one is much smaller than the other, but it’s not entirely clear to me how valid that would actually be at the CMB.
To clarify what you have said, basically step 5 is carried out first followed by (1+4), which is to say it is possible to insert steps 2 and 3 in between them? And I forgot to mention another reason why the Stokes’ equation is being solved twice is also because I need to compare u and v using some postprocessing function at each time step (e.g. comparing the relative magnitudes etc). Hence, I initially thought this was a more direct way to obtain the separate velocity fields for that purpose.
Kang: No, ASPECT doesn’t separate the solves for the velocity field that results from different effects. We just compute all forces that contribute to the flow field and only solve once. That’s the workflow I suggest to stick with because that’s what ASPECT already does. In your case, that might mean to compute boundary values u_{z,0} first and use this in the definition of boundary values for the one, global u+v/P_u+P_v. I recognize that this doesn’t allow you to later compare u and v because you only compute one velocity field.
Timo outlined some ideas of how you could add a second solve step, but ASPECT really isn’t set up to do that right now. It would require a substantial development effort. If you’re interested in this kind of thing, I’d take a look at the classes that compute melt generation/transport and the free surface module. Both of these extend ASPECT by additional equations, and that’s what you want to ultimately do. Neither of these were trivial additions to the ASPECT code base, and I would expect what you want to do to fall in a similar category.