geodynamics / aspect

A parallel, extensible finite element code to simulate convection in both 2D and 3D models.
https://aspect.geodynamics.org/
Other
223 stars 235 forks source link

Changing particle logic for iterated advection schemes? #4123

Closed anne-glerum closed 3 years ago

anne-glerum commented 3 years ago

There are potentially 2 locations where I think the logic regarding particles should be updated now that they can also be used with iterative Advection schemes and are advected at the beginning of the timestep.

  1. In solver_schemes.cc assemble_and_solve_composition: Do we need to advance the timestep at all, or only interpolate to the composition fields during the first timestep? The behavior now would be different for single and iterated Advection schemes I think. In the first case, velocity would be zero anyway (it's before the Stokes solve), and advection wouldn't change the particle positions. For iterative Advection schemes, velocity is not necessarily zero in subsequent iterations, so the particles could move. https://github.com/geodynamics/aspect/blob/7fbbe8083d9bc95f0f9fab6f32018f53c8b31311/source/simulator/solver_schemes.cc#L200-L201
  2. In initial_conditions.cc interpolate_particleproperties: Should the setting of old and old_old_solution only happen during the first nonlinear iteration? Should the current_linearization_point also be set, like for the compositional fields? https://github.com/geodynamics/aspect/blob/7fbbe8083d9bc95f0f9fab6f32018f53c8b31311/source/simulator/initial_conditions.cc#L337-L341
anne-glerum commented 3 years ago

After discussion with Rene:

  1. The particles are advected with a time step length of 0 during the zeroth time step during advance_timestep, so this advection has no result. However, advance_timestep also includes the updating of the particle properties, which is a necessary step. So no change needed.
  2. old_solution and old_old_solution should indeed only be set during the first nonlinear iteration at each initial mesh refinement level. I'll make a PR for that. The current linearization point is set by the calling function afterwards after all fields are interpolated, so no need to do that here.
gassmoeller commented 3 years ago

Closed by #4167