Support 180 degree phi-periodic 2d spherical shells

[Francyrad]

Summary

• Add support for phi-periodic 2d spherical shell geometries with a 180 degree opening angle.
• Keep 3d 180 degree spherical shells out of this PR, because that case needs separate benchmarks and documentation for periodicity, flow constraints, and nullspace handling.
• Update the Phi periodic parameter documentation to list only the supported 2d opening angles.
• Add regression tests for 2d periodic half shells and particle advection across the 180 degree periodic boundary.

This addresses the 2d spherical-shell part of #6926.

Testing

• git diff --check
• cmake --build build-debug --target aspect -j20
• cmake --build build-debug --target setup_tests -j20
• ctest -N -R 'periodic_half_shell_3d' -> Total Tests: 0
• ctest -R '^(periodic_half_shell|particle_periodic_half_shell|particle_periodic_half_shell_2|periodic_quarter_shell|particle_periodic_quarter_shell|particle_periodic_quarter_shell_2)$' --output-on-failure
• Result: 100% tests passed, 0 tests failed out of 6 with numdiff enabled for numerical comparisons.

[gassmoeller]

I know this is just a draft for now, I just wanted to mention a few concerns I have about this PR.

The periodicity of a 3D spherical shell is fundamentally different from what we have implemented at the moment for 2D. In a 2D spherical shell the periodic boundaries are not physically at the same location, there is no location where the periodic faces that are connected across the periodic boundary connect to each other. In a 3D 180 degree spherical shell this is different, the two periodic faces (the two half cuts across the sphere that form a plane) are essentially just one face that you split into two periodic faces. This has a few consequences:

1. It is not unique how to make this face periodic. You chose the phi direction here (at least from the description of the PR), and this is ok, you have to choose something, but this needs to be documented for the user.
2. How you choose the periodicity direction matters, because it will have consequences for the movement of material. Your periodicity matrix says for example that the velocity normal to the periodic face needs to be equal to minus the velocity on the other periodic face. If you now think about the vertex where the two faces touch, this can only be true if the velocity is constrained to 0. So I assume your periodic boundary really implies a no velocity constraint along an axis through your model. This is not obvious for someone activating this boundary. And it is also not clear if deal.IIs periodicity constraints will actually enforce this (or something else that is wrong entirely).
3. Adding a 180 degree 3D shell even without periodicity adds a rotational nullspace to your model solution. The model could rotate without resistance around a spin axis from the origin to the center of the shell. It is not clear our nullspace rotation handles this correctly at the moment.

I think what I am saying is mostly, this could be an interesting contribution, but it will need:

• quantitative benchmarks that show the new functionality applies 'correct' periodic boundary conditions and which consequences this has for the flow. This can for example be achieved by prescribing known flow fields (like solid body rotation) in parts of the domain, and checking that these flow field correctly transport across the periodic boundaries
• benchmarks that show nullspace removal works correctly for the 180 degree spherical shell with and without periodic boundary conditions
• a cookbook with documentation on what a periodic boundary means for a 180 degree spherical shell and which assumptions this requires

[Francyrad]

Adding only 180° in 2D for the moment

[Francyrad]

I updated this PR to keep only the 2D 180 degree phi-periodic spherical shell support.

The 3D 180 degree case has been removed from the scope of this PR, because I agree that it needs separate benchmarks, documentation, and checks for the periodicity direction and rotational nullspace behavior.

So the current PR is intentionally limited to the simpler 2D case, analogous to the existing 90 degree periodic shell tests, with added 180 degree periodic shell and particle advection tests.

Would this reduced 2D-only scope be acceptable for review/merge, while leaving the 3D 180 degree case for a separate future contribution?
@gassmoeller

[gassmoeller]

This should work in principle, but I have some comments that I would like to see addressed before merging this.

Also: Thanks for adding the test. Did you check if particles are crossing the periodic boundary correctly, i.e. did you see a particle with a certain ID leave the domain on one side and enter it on the other? I am just asking because from the setup I cannot tell if this actually happens in your test model.

[gassmoeller]

I dont think this function adds anything to the code for the 2D version (maybe there was a use in 3D). Please remove.

[gassmoeller]

Like I mentioned above, I think the previous code was clearer here. Please revert this line.

[gassmoeller]

I dont understand this piece of code. It looks like you are rotating the whole grid by 90 degrees. Is this necessary or a preference? It should definitely get a comment explaining the reasoning behind it. But this also creates practical problems. Consider a model that you first run and evaluate without periodic boundaries. Then you switch the boundary condition and suddenly your model geometry changes. This is unexpected from a user perspective and may break workflows that assume they can rely on the geometry of the 180 degree 2D spherical shell. Does the periodic boundary work without rotating the grid?

[gassmoeller]

So this part relies on the geometry modification you did above. But couldnt this be formulated no matter how the geometry is rotated?

Also I would recommend to keep the structure in a way that there is only one early return statement, i.e. the meaning of this condition should be:

  if (phi == 180 && coordinate suggests the position needs to be rotated)
    set rotation matrix
  else if (phi == 90 && coordinate suggests the position needs to be rotated)
    set rotation matrix
  else
    // nothing to do, return early

  // now do the rotation

[gassmoeller]

This code now looks odd, because you introduced the function phi_periodicity_rotation_tensor above, which should already give the rotation matrix. Instead you could write the following to make the intent more clear:

Suggested change

	else if (position[0] < 0.)
	{
	// if the position crossed the west boundary there
	// is nothing to do, use the rotation matrix as is
	// to rotate the position clock-wise
	}
	else if (position[1] < 0.)
	{
	// rotate the opposite direction (anti clock-wise)
	// if the position crossed the east boundary
	rotation_matrix *= -1;
	}
	else

[gassmoeller]

It is unfortunate that GridTools::collect_periodic_faces requires a FullMatrix as input, while we want to use a Tensor<2,dim elsewhere, but the solution is not to duplicate this code. Please check if there is a deal.II function to convert a dimxdim matrix into a Tensor<2,dim>, and if not, then we can easily write our own:

template <int dim>
FullMatrix<double>
phi_periodicity_rotation_matrix(const double phi)
{
  FullMatrix<double> rotation_matrix(dim);
  const Tensor<2,dim> rotation_tensor = phi_periodicity_rotation_tensor(phi);

  // loop to convert rotation tensor into rotation matrix
  ...

This avoids the code duplication.

[gassmoeller]

Also please add a changelog entry.

[Francyrad]

This should work in principle, but I have some comments that I would like to see addressed before merging this.

Also: Thanks for adding the test. Did you check if particles are crossing the periodic boundary correctly, i.e. did you see a particle with a certain ID leave the domain on one side and enter it on the other? I am just asking because from the setup I cannot tell if this actually happens in your test model.

Hi, for my simulations cases (LLSVPs simulations, so thermochemical simulations where i give a field for LLSVPs and a field for the mantle) LLSVPs cross the boundary perfectly like in the 90° cases without any problem, I was also thinking that these 180° implementation was "good enough" to "replicate" the behaviour of the 360° sims. I'm gonna address also your comment and make clearer tests

[Francyrad]

Hi @gassmoeller, I went through your comments and cleaned up the patch.

What I changed:

• I removed the extra helper that only returned the periodic direction.
• I removed the extra rotation of the mesh for the 180 degree case. This means that switching Phi periodic on does not silently rotate the geometry anymore.
• I kept one canonical rotation tensor and only convert it to FullMatrix where deal.II needs that type.
• I rewrote the particle periodicity handling so the 180 degree case is separate from the 90 degree case.
• I updated the affected particle test references.
• I added a changelog entry.

The relevant local tests for 2d spherical-shell periodicity and particles pass.

I also did a larger manual sanity test with a Steinberger-style 180 degree shell and particles. I first ran it without phi periodicity to 100 Myr. The 100 Myr non-periodic output had 43,425 particles. Then I restarted from that run with Phi periodic = true; the restarted case produced outputs through 1 Ga without crashing. The first particle output after the periodic restart had 43,585 particles, and the final 1 Ga output had 63,862 particles.

The particle count is not expected to stay constant in this setup, because the model uses the remove and add particles load-balancing strategy with minimum and maximum particle limits per cell. Therefore ASPECT can add particles during the run when cells become under-populated; the important check here is that particle output continues correctly after restart and that the particle fields are present.

To make the restart behavior easier to check, I also made a side-by-side movie from the original output files. The left panel shows the particles colored by temperature, the right panel shows the full temperature field, and the frame around 100 Myr is marked with Restart: Phi periodicity enabled.

The zipped movie is available here:
https://github.com/Francyrad/aspect_dynamic_core_fixes/raw/pr-6957-video-assets/pr-assets/6957/steinberger_180_100M_nonperiodic_then_1Ga_periodic_english_slow.zip

[gassmoeller]

Very nice, just two minor comments left. Thanks for checking, I just couldnt tell from your previous comments if you had actually observed particles crossing the boundary.

After addressing the comments and updating the test results this should be ready.

[gassmoeller]

Sorry, this file is reserved for changes that happened before 3.0.0. Please put this changelog entry into a new file in /doc/modules/changes. We will combine all changelog entries into a new file for the release 3.1 later. (see here for instructions).

[gassmoeller]

The two particle tests run with a global resolution of 4, which makes the tests run in less than 10 seconds, which saves some compute time. Since you have to update the test results anyway (the tester failed), can you also reduce the resolution here to 4?

instructions on updating test results: https://aspect-documentation.readthedocs.io/en/latest/user/extending/testing/writing-tests.html#updating-test-results-for-pull-requests

[Francyrad]

@gassmoeller done

[gassmoeller]

Thanks, then this is ready when the testers are done.

[gassmoeller]

/rebuild

[gassmoeller]

One of our testers has connection issues, but it is almost a duplicate of one of the other test setups, so I think this is fine to merge now.