Reorder particle advection for use with mesh deformation

[anne-glerum]

This commit

1. removes the Assert that disables using particles icw mesh deformation
2. introduces a post_mesh_deformation signal that the particle_world connects to to sort the particles after mesh deformation;
3. advects particles before assemblying and solving composition, instead of at the end of the timestep;
4. adds a test that checks the reference position of a single particle;
5. adds a particle property that outputs the reference position of the particle.

DONE: check particle advection is correct for the iterative solver schemes.
DONE: add tests
DONE: change test that checks that particles + mesh deformation fails
TODO: update test results

For all pull requests:

• I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

• I have tested my new feature locally to ensure it is correct.
• I have created a testcase for the new feature/benchmark in the tests/ directory.
• I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

[gassmoeller]

I think this looks very good! The only major concern I have is for the particle properties. Most properties are fine (everything that does not actually store history, but only extracts information at the current state does not care), but the ...strain... properties and the pT path property will probably be wrong in the iterated schemes. Could you try storing a copy of the ParticleHandler and resetting that?

[gassmoeller]

Add another line (3 empty lines is our default between functions in .cc files, 1 line in .h files).

[anne-glerum]

So what to do when there is a mix of 1, 2, and 3 empty lines in the .cc file (as in this case)? I shouldn't adjust them all in this PR right?

[anne-glerum]

I've made the iterated Advection solver schemes still use the old order of particle advection at the end of the timestep. When we have the functionality to reset the ParticleHandler, we can change those 2 schemes as well.

[anne-glerum]

Hi @gassmoeller , I've used your copy_particle_handler function (#3818) to also reorder the particle advection in the iterated Advection schemes, thanks! I tested the schemes with the test we set up with @ricitron and with a similar setup that does not deform the mesh but prescribes a non-zero Stokes velocity. In both cases particle positions and reference positions are correct. A van Keken benchmark also runs with a larger number of particles.

Could you have another look? Also, a lot of tests will fail due to the reordering of the particle advection. In how far do we have to check whether the diffs make sense? Are there tests with known solutions to check?

[gassmoeller]

Just some minor comments, in general this looks good. In terms of benchmarking: I think the best benchmark will be benchmarks/rigid_shear and the simplified test case rigid_shear_time_dependent. This is the time dependent box benchmark from our paper (Gassmöller, R., Lokavarapu, H., Bangerth, W., & Puckett, E. G. (2019). Evaluating the accuracy of hybrid finite element/particle-in-cell methods for modelling incompressible Stokes flow. Geophysical Journal International, 219(3), 1915-1938.)

From the test results it looks like the output changes only slightly (because we now use the extrapolated velocity instead of the real velocity?), but it could be good to confirm by running the full benchmark again (cheap) and not only the test, and we should also run the same benchmark with a nonlinear solver scheme. I am working on getting another benchmark into a PR, but the rigid shear one should be sufficient already.

[anne-glerum]

From the test results it looks like the output changes only slightly (because we now use the extrapolated velocity instead of the real velocity?), but it could be good to confirm by running the full benchmark again (cheap) and not only the test, and we should also run the same benchmark with a nonlinear solver scheme. I am working on getting another benchmark into a PR, but the rigid shear one should be sufficient already.

@gassmoeller. I've now run part of the rigid shear time dependent benchmark (lowest and highest resolution, only for cell average interpolator and for 10, 45 or 80 particles per dim). The change in errors is small and decreases with resolution and with the number of particles:

I'll run a different suite with a nonlinear solver scheme.

[anne-glerum]

Hi @gassmoeller . Here results for the benchmark using the iterated Advection and Stokes scheme and bilinear least squares. The differences are small; but biggest for Q3Q2. I didn't run the complete suite of models, should I?

[gassmoeller]

The results looks great, thanks! No need for further tests from my point of view, if something would be wrong it should show up prominently in these plots. Can you update the test output and add a changelog entry? Afterwards this should be good to go.

[gassmoeller]

Sorry to bring up one more point: Could you check what is happening in the viscoelastic tests (e.g. visco_plastic_vep_stress_build-up_yield_particles)? It seems like the properties are now updated one timestep earlier than before, i.e. compositions were 0/0/0 before, and now 1e6/1e6/.... The values are the same as in later timesteps, it just looks like the update happens one timestep earlier than before. Maybe that is even correct, but I would like to understand why it happens.

[anne-glerum]

Hi @gassmoeller good point! I haven't quite figured it out yet. This is just to help myself remember:

• These tests store the stress tensor elements on compositional fields, which are changed through reaction terms set in the material model, and advected with the particle method.

• In the initial timestep, velocity is assumed zero before solving the Stokes equation.

• In the initial timestep, reaction outputs are set to zero by the function fill_reaction_outputs.

• On mainline, in the initial timestep, particles are advected at the end of the timestep, and the compositional values have not been changed from zero (advection equation not solved + velocity = 0 + reaction term = 0).

• In this PR, in the initial timestep, particles are advected before the advection equation with a velocity of zero, and the reaction terms are also zero.

• old_solution and old_old_solution are set to the t0 solution

• On mainline, at t1, we solve for velocity and then advect the particles and update the elastic stresses through the reaction terms, which depend on the current solution (that of t1).

• In this PR, at t1, we advect the particles, update the elastic stresses through the reaction terms which depend on the current solution (that of t0), and then solve for velocity.

[anne-glerum]

This is what happens for the failing vep tests:
At t0, reaction_terms are not evaluated, so the elastic stresses stored on the particles remain zero.
At t1,

1. Particles are advected and updated, which includes ElasticStress::update_particle_property
2. Particle properties (elastic stresses) are interpolated to compositional fields → old_stress no longer zero
3. Stokes assembly using the compositional fields → also filling the stress force rhs terms and the reaction terms based on the old_stress, which is the current composition of t1
4. Solve Stokes with the elastic force rhs term
5. During postprocessing, Elasticity<dim>::fill_elastic_force_outputs and Elasticity<dim>::fill_reaction_outputs are again evaluated, but not needed. Setting in.requests_property(MaterialProperties::reaction_terms) could prevent that?

First I thought old_stress should be the stress from the last timestep, but without particles, the old_stress would also be the compositional value of the current timestep. So I think it's correct, but I'd like to discuss this with @gassmoeller and @naliboff . I noticed the bending beam benchmark has a field and a particles prm, I'll compare the two.

Previously, the particles were advected at the end of the timestep. So at t0, the update of the property elastic stress does nothing, as reaction terms are not evaluated. At t1, this zero stress property is interpolated onto the compositional fields, which are therefore also zero. At the end of timestep t1, the particles are advected and updated and carry nonzero stresses, while the stress compositional fields are still zero (see figure below showing t1=1000yr from the test visco_plastic_vep_stress_build-up_yield_particles).

[anne-glerum]

There are also a lot of tests where the particle velocity outputted in the gnuplot file for t0 is now 0. This is because the particles are now advected before there is an actual velocity solution, instead of at the end of the timestep.

[gassmoeller]

Thanks for checking, all of this looks correct and like an improvement over the current state. If you update the remaining tests I am fine with merging this.

During postprocessing, Elasticity::fill_elastic_force_outputs and Elasticity::fill_reaction_outputs are again evaluated, but not needed. Setting in.requests_property(MaterialProperties::reaction_terms) could prevent that?

Only if the material model also checks requests_property before calculating these properties. Currently this structure is more of a placeholder that can be used, but is not used in most material models (the default is to compute all properties).

[anne-glerum]

Hi @gassmoeller I've stored the particle_handler copy in world.cc, as requested in #3932. If this is what you were looking for, I'll update the test results.

During postprocessing, Elasticity::fill_elastic_force_outputs and Elasticity::fill_reaction_outputs are again evaluated, but not needed. Setting in.requests_property(MaterialProperties::reaction_terms) could prevent that?

Only if the material model also checks requests_property before calculating these properties. Currently this structure is more of a placeholder that can be used, but is not used in most material models (the default is to compute all properties).

The functions in elasticity.cc do check whether the reaction_terms are requested.

[naliboff]

@anne-glerum, @gassmoeller - the updated elasticity results look correct. In fact, this is an improvement. Due to the order of how the compositional field values for elasticity were updated when tracking elasticity on particles the output compositional field values always appeared to be a one time step off. These listed compositional field values between tests using fields or particles to track the elastic stresses are now in line. Thanks for your work on this and happy to discuss or test in more detail if needed.

[gassmoeller]

Yes, this is exactly what I had in mind. Go ahead and update the tests.

The functions in elasticity.cc do check whether the reaction_terms are requested.

Oh, yes, then not requesting reaction terms should stop the computation. Feel free to adjust.

[gassmoeller]

one empty line too much :-)

[anne-glerum]

@gassmoeller While updating the tests, I realized that the tests I added in this PR don't have the right results yet. Unfortunately the mesh_deformation_particles test doesn't work correctly anymore in the sense that the top boundary is not displaced through the function plugin. I don't think it's a problem with this PR, as it's the same with a master build I had lying around. The test has single Advection, no Stokes scheme, prescribed Stokes solution of zero, plus a top mesh boundary with a function prescribing upward motion of the surface. Will have to look into why the mesh is not displaced.

[anne-glerum]

Alright, fixed the test. However, mesh_deformation_particles currently does not produce the correct reference location. That's not because the particle reference location is incorrect after mesh deformation at the beginning of the timestep (we do a sort), but because the test post_mesh_deformation signal is evaluated before the solve_timestep post_mesh_deformation signal that does the sort. @gassmoeller What determines the order in which the signals are handled?

[gassmoeller]

Slots (the functions connected to the signals) are executed in the order they were connected. Unfortunately the test connects through the signal connector function (which is called from core.cc:149), while the particle world only connects later in core.cc:418. I see two ways to fix this, both not quite clean:

1. You could modify the test to have another function connect_to_signal that is connected by the signal_connector function to a signal that is triggered at a later time (e.g. the pre_set_initial_state signal). This way at the beginning only connect_to_signal is connected to pre_set_initial_state , and only when pre_set_initial_state is triggered is the test function post_mesh_deformation connected to the post_mesh_deformation signal (at which point the particle world has connected its function already).
2. You can change the connection of the particle world to enforce being the first slot that is called. This can be done by changing the connect call in world.cc:276 to (i hope this works, have not tested, only got this from the boost documentation):

     signals.post_mesh_deformation.connect(
        [&] (const SimulatorAccess<dim> &)
      {
        particle_handler->sort_particles_into_subdomains_and_cells();
      },
boost::signals2::at_front);

2. is the smaller change, and may be a good idea anyway to ensure the particle sorting is always the first thing happening after mesh deformation.

[anne-glerum]

Option 2 worked like a charm, thanks @gassmoeller! Waiting for the updated test results and then this should be done :)

BTW, I had to patch the indentation with the tester results, as my local ninja indent didn't update what the tester thinks is wrong. Does that happen more often?

[anne-glerum]

@gassmoeller Tests are updated, should be ready.

[gassmoeller]

Looks like one of the tests is still failing.

[anne-glerum]

This test now passes on my machine. Are the differences with the tester ok to just update the results?

[gassmoeller]

Yes, if the test works for you then just update the results from the tester.

[gassmoeller]

Very nice, thanks for the changes!