synchronize particle movement and compositional fields

[gassmoeller]

Currently the composition and temperature advection equations are solved before the Stokes solution, but the particle advection happens afterwards in the post processing stage. I would prefer that all advection processes happen at the same time, this prevents surprises in the case of active particles. For passive particles it does not matter at all.

[gassmoeller]

See #828 for some discussion about this issue.

[bangerth]

But they already don't happen at the same time. We first solve the temperature equation, then solve the composition equation for composition 1 with the now current temperature but the old values for all other compositions, etc. I think this sort of thing is unavoidable in an operator splitting approach unless we want to iterate things out. So I don't see this as a significant problem in search of a solution. Running all of this as a postprocessor also has the advantage of being simple.

But that's just one opinion. What do others think?

[gassmoeller]

I see your point of why it is simpler to keep everything in the postprocess stage. However, the argument that the temperature and composition are already updated right after their individual solution does not quite convince me. Usually there is little interaction between the solution of these advection equations, the important part is that the Stokes system gets an updated estimate for the density. So if you imagine an interface that in some sort represents a density jump, you can currently realize this by temperature field, composition field or particles equally well. However the velocity solution for the two field methods will be identical, while the solution with particles would be "off" by one timestep. I think this will cause a lot of trouble for myself (comparing models will become a nightmare) as well as for other users, who will simply attribute this to the difference in advection method, while in reality it is a difference in the order of our time-stepping. But I would like to hear other opinions. Maybe @tjhei or @jdannberg have suggestions?

[jdannberg]

Okay, here is my opinion: 1) Modellers in geodynamics have been using particles for a long time, and that is what they are used to. We now have both fields and particles, so people will want to compare. So we should also implement these features in a way that users can easily compare between the two approaches and see which one works better for their applications. 2) Once we have active particles, I would find it confusing if they were still postprocessors. They will actively change the behaviour of the model, so they should probably become something else than a postprocessor. 3) Particles potentially advect similar properties as compositional fields, and could also be used to advect grain size/porosity etc. So at some point we may want to iterate things out, so it would be useful to have them somewhere inside of the iterations, and not in the postprocessing step. 4) I also see the point Rene makes: the temperature usually does not affect the advection of the compositional fields/particles, and the compositional fields might influence each other, but they are only updated after all of them are solved. So I think it would make sense to advect all quantities we want to advect with the same solution for the velocity.

[tjhei]

2) Once we have active particles, I would find it confusing if they were still postprocessors. They will actively change the behaviour of the model, so they should probably become something else than a postprocessor.

Agreed. Compositional fields are also not postprocessors. :-)

Rene, you would advect particles after temperature and compositional fields, but before Stokes? Then you no longer have access to the current Stokes solution?!

[egpuckett]

I have just read this discussion, meaning I haven't spent hours or days thinking about it. However, my first thought is that it is more consistent with the structure of the overall model to advect active particles at the same time one advects the temperature and compositional fields. One line of reasoning is that the temperature and (possibly) the compositional field(s) will change the force on the RHS of the Stokes equations and also possibly the viscosity, etc. when one advances the velocity from time t{n} to time t{n+1}. If one were to do an error / accuracy analysis, it seems to me that if the particles also contribute to a change in the Stokes equations that are being solved to advance the velocity from time t{n} to time t{n+1}, then it would be easier to do the analysis if the change due to the particles occurs at the same time as the changes to the Stokes equations due to the temperature and compositional field(s). In fact, I can imagine such an error analysis not working out if there were a difference in times for temperature / composition versus active particles.

Another thought is that it might make sense to review some of the literature on 'particle in cell' methods, especially articles that include analysis or convergence studies.

Perhaps it would be helpful to write the equations out for a simple case and study them. Is there a simple case of active particles in the geophysical literature that we could examine / think about?

Finally, I think @jdannberg's argument:

" We now have both fields and particles, so people will want to compare. So we should also
implement these features in a way that users can easily compare between the two approaches and see which one works better for their applications."

makes a lot of sense. If the active particles' properties are applied in the post processing stage, but the compositional field's properties are applied before the Stokes solve, one could not compare the efficacy of using particles carrying property A against using a compositional field carrying property A. It would be like comparing apples and oranges, as they say.

[bangerth]

I see - I will be happy to submit to the wisdom of the many! :-)

[egpuckett]

This is truly what I treasure about the ASPECT project. There is an egalitarian approach in which each participant's view is valued, if not viewed equally. I've worked on similar projects where the views of the project leader were sacrosanct; i.e., "My way or the highway". This is what sets the ASPECT project apart from many similar efforts.

[bangerth]

Thank you. I appreciate the kind words!

[egpuckett]

My pleasure!

[egpuckett]

I occurs to me that from the point of view of the routines euler.cc, rk_2.cc, and rk_4.cc, if one simply calls them before the Stokes solve, the velocity is what is currently old_velocity; namely the velocity at time t_n, then the particle movement and compositional fields will be updated at the same point in the computation. However, I imagine that there is a great deal more to turning a postprocessor into a part of the main computation. Am I missing anything?

[egpuckett]

Having giving the matter some thought, if we advect the particles before the Stokes solve, we cannot use rk_2.cc, and rk_4.cc as particle integration methods (i.e., to update the positions of the particles in time). Rene, what were you planning to use as an ODE solver to advect the particles before the Stokes solve?

[bangerth]

@egpuckett -- you won't be able to do it in the first time step, but you certainly have all of the information about the previous time steps after that. Or do I misunderstand?

[egpuckett]

I don't think we'll have the information after the first time step for the following reason. If we move the particle integration step to before the Stokes solve, then (velocity) will be the velocity at time tn, not t{n+1}. Our current implementation of, say, rk_2.cc relies on having (velocity) be the velocity at time t{n+1}. The reason we need the velocity at t{n+1} is that we approximate the velocity at time t_{n+1/2} by forming the average

( (velocity) + (old_velocity) ) / 2 ,

which at postprocessor time is an approximation to

( v(tn, x) + v(t{n+1}, x) ) / 2 v(t_{n+1/2}, x)

We can't do that before the Stokes solve.

We could use BDF2, or an explicit multistep method. In the work @yinghe616 has done (see PR #864) she has used an explicit Strong Stability Preserving (SSP) method; e.g., see

Gottlieb, Shu, Tadmor (2001) SIAM REVIEW Vol. 43, No. 1,pp. 89–112

outside of ASPECT, since apparently the BP limiter is more effective with an explicit ODE solver. However, I believe the implementation in PR #864 uses BDF2 and that in a future PR she will implement the SSP ODE solver.

By the way, I'm telling everyone here at UCD / CIG including @hlokavarapu, @yinghe616, and @naliboff what you wrote to me recently:

Rule #1: Only one functional change per pull request or commit! :-)

So, if you suddenly see a pile of PRs from us, that is the reason. Apparently, most of us have the tendency not to put our code on display until it is 'perfect' and 'complete', or something along those lines.

[bangerth]

In the first time step, we solve the equations with dt=0 anyway (we solve for the advection equations first, but we don't have a velocity yet, so dt=0 is the only reasonable choice).

But in general, I see no reason why we shouldn't use something like an explicit Euler scheme for the first time step, and only BDF2/RK3/RK4 for successive time steps. This is what we do for the advection solvers already, IIRC.

[yinghe616]

@egpuckett , I think @bangerth mentioned above is that at the first time step 0, aspect doesn't compute the velocity field. In the codes, actually the function call compute_initial_pressure_field (), (after set_initial_temperature_and_compositional_fields () ), will set all the velocity fields all zero values; And then it moves to the next time step 1, which is different than we thought before - solving stokes first in order to get the initial velocity at t=0 after initializing the temperature and composition.

Below is the code at time step 0: set_initial_temperature_and_compositional_fields (); compute_initial_pressure_field (); initialize_tracers ();

Therefore, at any time step k, the computation order will be : solve T -->solve composition C -->solve stokes u+p and move to next time step k+1

So, any post-processor calls should have the numerical solutions T, C, u, p at the same time.

[bangerth]

Almost correct. In the first time step, we don't solve for T, C (or, better, we do solve, but with zero time step and therefore remain at the initial conditions), but we definitely solve for the Stokes solution (which we need to do, since no initial conditions are given for u/p, so we need to do a solve to find them everywhere).

This is correct:

Therefore, at any time step k, the computation order will be : solve T -->solve composition C -->solve stokes u+p and move to next time step k+1

So, any post-processor calls should have the numerical solutions T, C, u, p at the same time.

[yinghe616]

@bangerth, Thanks for the clarification. It makes more sense now. :p

[egpuckett]

But in general, I see no reason why we shouldn't use something like an explicit Euler scheme for the first time step, and only BDF2/RK3/RK4 for successive time steps. This is what we do for the advection solvers already, IIRC.

Do you have and explicit RK2 and / or RK3 & RK4 currently implemented for the advection of a compositional field? I don't see you one can do that, since one needs the value of the velocity at the half time level v(t_{n+1/2}, x) for, say, RK2, which can only be obtained by extrapolation if the velocity is only known at

(velocity) = v(t_n, x), (oldvelocity) = v(t{n-1}, x) ) and (*old_oldvelocity) = v(t{n-2}, x)

[yinghe616]

@egpuckett, current aspect solves the advection term implicitly (or semi-implicitly: term adaptive BDF2 for dC/dt + u^k\cdot \nabla C^{k+1}, for time step k+1, where u^k is from the previous step k, and so it's not necessary to do the extrapolation as long as C uses C^{k+1}). @bangerth, please correct me if I was wrong.

[yinghe616]

Just realized your concern that u{k-2} will be overwritten by u{k-1} once we have u{k+1} as we can only store three solutions. Then why don't we just use the average (u{k+1}+u{k})/2 to approximate u{k+1/2}? Is there anything I am still missing?

[egpuckett]

You one can't do that, since one needs the value of the velocity at the half time level v(t_{n+1/2}, x) for, say, RK2, which can only be obtained by extrapolation if the velocity is only known at

(_velocity) = v(t_n, x), (_oldvelocity) = v(t{n-1}, x) ) and (*old_oldvelocity) = v(t{n-2}, x)

This is the case when the compositional & tracer advection is done BEFORE THE STOKES FLOW.

On 5/10/2016 12:06 PM, Ying He wrote:

Just realized your concern that u{k-2} will be overwritten by u{k-1} once we have u{k+1} as we can only store three solutions. Then why don't we just use the average (u{k+1}+u{k})/2 to approximate u{k+1/2}? Is there anything I am still missing?

— You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub https://github.com/geodynamics/aspect/issues/829#issuecomment-218258134