Fix mesh deformation for GMG by adding support for higher-order FE degrees

[tiannh7]

Hi all,

This PR fixes mesh deformation issues in GMG by enabling support for higher-order finite element degrees, addressing anomalies at mesh refinement boundaries discussed in issues #6835 and #6248.
It extends the improvements introduced in PR #6284 related to free surface accuracy and mesh deformation.

[bobmyhill]

Great, thanks! This will need someone else to take a look over the changes, but I've left some basic comments above.

A few other questions:

• Is this a standalone PR, or does it need #6284 to be merged?
• Could you create a test, or point us to an existing test that fails on main and is fixed here?
• Can you create a changelog for this fix?

[bobmyhill]

Can you use whole words in your function names, and make them descriptive?
Doing this will help speed up code comprehension.

[bobmyhill]

Should this be a user parameter, or is there always a "correct" answer for a given problem?

[tjhei]

For backwards compatibility, we probably want the default to be 1. Using the velocity degree (so 2 for most computations) is probably the best approach. That would mean we need an option

[bobmyhill]

If this should be a user parameter, can you add more documentation here? Remember, this is probably all the documentation that most people will see, so it needs to be sufficient to give users an idea of the values they should use (or another place to look for the information).

[tiannh7]

Hi @bobmyhill,

Thanks a lot for your valuable comments and questions — they are all very helpful and I completely agree with your points.

Is this a standalone PR, or does it need #6284 to be merged?

This is a standalone PR. Many thanks to @gassmoeller for the foundational work in #6284!

I plan to address the review comments, add the requested test, and prepare a changelog by tomorrow afternoon or at the latest within this week.

[tjhei]

Suggested change

	throw std::runtime_error("Unsupported mesh deformation polynomial degree!");
	AssertThrow(false, ExcMessage("Unsupported mesh deformation polynomial degree!"));

[tjhei]

For backwards compatibility, we probably want the default to be 1. Using the velocity degree (so 2 for most computations) is probably the best approach. That would mean we need an option

[tjhei]

Awesome, thank you for working on this!

[tiannh7]

Thank you both very much for your detailed reviews and valuable feedback! I have addressed the requested changes.

[bobmyhill]

The tester failures are all expected.

Could you add the new parameter to all but one of the failing tests, and leave boundary_traction_function_cartesian_free_surface as a test of the assert (i.e., just patching the test output for that failed test)?

[tjhei]

I don't quite understand this. Why do you need "Use stokes velocity polynomial degree" if the user can just specify the value?

[tiannh7]

I will remove the redundant option Use stokes velocity polynomial degree for mesh deformation.

[tiannh7]

Hi all,

I ran a small automation script to batch-run the test suite and refresh reference outputs. As a result you may see a number of files added, modified, or removed in this PR (these are mostly regenerated test outputs). If anything here looks incorrect or unexpected, please let me know which files to revert or investigate and I’ll address it promptly.

Note: while running the tests I observed a non-convergence in the geoid_freesurface case — see the excerpt below:

----------------------------------------------------
Exception 'ExcMessage (exception_message.str())' on rank 0 on processing:

--------------------------------------------------------
An error occurred in line <2844> of file utilities.cc in function
    void aspect::Utilities::throw_linear_solver_failure_exception(const std::string &, const std::string &, const std::vector<SolverControl> &, const std::exception &, const MPI_Comm, const std::string &)
Additional information:
    The iterative Stokes solver in
    StokesMatrixFreeHandlerLocalSmoothingImplementation::solve did not
    converge.

    The initial residual was: 4.365282e-02
    The final residual is: 4.965443e-07
    The required residual for convergence is: 4.365282e-09
    See output-geoid_freesurface/solver_history.txt for the full
    convergence history.

    The solver reported the following error:
----------------------------------------------------

I have previously seen this behave like a solver issue. Once this PR is essentially complete I will open a dedicated follow-up thread to propose a fix.

[bobmyhill]

To update test results, you can use the patch .diff from the 9.6 tester (there's a link in the tester output, e.g. here https://github.com/geodynamics/aspect/actions/runs/21020930699/artifacts/5136857595). You can apply the patch with git apply <filename>

[tiannh7]

@bobmyhill Thank you!

[tjhei]

Suggested change

	this->get_pcout() << " Solving mesh displacement system... ";
	this->get_pcout() << " Solving mesh displacement system... " << std::flush;

Otherwise, the output will likely not be printed until you hit the std::endl

[tjhei]

This is worth thinking about:
By doing this check you make sure people get accurate results. On the other hand, you are breaking existing simulations.

[gassmoeller]

I am ok either way. The Assert seems slightly too cautious, but using a Stokes degree 1 is quite exotic anyway (since it needs additional stabilization for the finite element to even solve the Stokes equation). So everyone who knows what they are doing and select a velocity degree of 1 likely also knows how to comment out this assert if they are convinced that is what they want.

[tjhei]

same here

[tjhei]

You have to do this now because you are changing the mapping at a later time now, right? Can you leave a comment here about this please? @gassmoeller is this problematic/expensive?

[gassmoeller]

I discovered the need for this in #6284 and did not extensively test or measure it. But: It is done before any particles are created and at this point it is just associating the particle handler with this mapping and triangulation. So at worst we duplicate the initialization of one (or more) empty particle handlers. I wouldnt expect it to be expensive.

[tjhei]

and Timo Heister :-)

[tiannh7]

Many thanks for the valuable review!
Your input has significantly improved the quality of this PR.

[tjhei]

I thought about this some more and think we could also decide to always use the Stokes degree:

1. I expect the increase in computational cost to be small.
2. It is more accurate and every user will automatically get this advantage instead of having to know that they need to update (or getting an error prompting them to do the change).
3. It is simpler.

If you want to allow the option to still use degree=1 I would suggest the following parameter and logic:
Mesh deformation polynomial degree has default 0 with the meaning of automatically picking an appropriate degree. I would use the Stokes degree in that case but you could also let 0 to mean "Stokes degree if curved, 1 otherwise".

What do you think?

[tjhei]

Note, test boundary_traction_function_cartesian_free_surface fails with:

An error occurred in line <241> of file </jenkins/workspace/aspect_PR-6836/source/mesh_deformation/interface.cc> in function
    void aspect::MeshDeformation::MeshDeformationHandler<dim>::initialize() [with int dim = 2]
The violated condition was: 
    this->get_parameters().mesh_deformation_polynomial_degree > 1

[tiannh7]

I thought about this some more and think we could also decide to always use the Stokes degree:

1. I expect the increase in computational cost to be small.

2. It is more accurate and every user will automatically get this advantage instead of having to know that they need to update (or getting an error prompting them to do the change).

3. It is simpler.

If you want to allow the option to still use degree=1 I would suggest the following parameter and logic:

Mesh deformation polynomial degree has default 0 with the meaning of automatically picking an appropriate degree. I would use the Stokes degree in that case but you could also let 0 to mean "Stokes degree if curved, 1 otherwise".

What do you think?

These are all excellent points! @bobmyhill @gassmoeller I’d love to get your take on which direction we should go.

[bobmyhill]

Hi @tiannh7,

I like the idea of always using the Stokes degree - this was the internal decision making I was thinking about in my first review, and it's great to learn that there is a reasonable choice. Thanks @tjhei!

FWIW I'm not a fan of "Stokes degree if curved, 1 otherwise" - it makes the logic more complicated and once the mesh is deformed, the distinction seems artificial to me.

I'm not strongly against keeping the new parameter, but I am also not strongly in favour if the only purpose is backwards compatibility; everyone's prm files should still work, they'll just work better.

[tiannh7]

Thank you both!
Changed mesh deformation to use the Stokes velocity degree by default (removed the separate parameter). Please report any issues. I’ll apply follow-up diffs after tests finish.

[gassmoeller]

Hi all, thanks for the improvements @tiannh7. I dont have time for detailed review today, I will take another look tomorrow. I am strongly in favor of the higher order mesh deformation FE (and I agree choosing the same degree as the stokes velocity makes sense). The reason I didnt progress further with my earlier try in #6284 was that I had an example model that didnt seem to be stable with the higher order surface deformation (developing oscillations after a while). However, maybe that is fixed here, or I didnt check correctly. I will see if I can find my old setup and test it against this branch. I will let you know when I have results.

[gassmoeller]

So I looked through the code and also ran some test models. I marked up some questions I had, they are mostly concerned with the mathematical background.

• 1. I think we need to identify good degrees for the mapping. Intuitively I would think this->get_geometry_model().has_curved_elements() ? std::max(4,mesh_deformation_fe.degree) : mesh_deformation_fe.degree would make sense. However this is different from the current choice, and I dont know how much testing you already did @tiannh7, maybe your choice works better? Why did you choose the numbers you have at the moment?
• 2. Then we should make sure the GMG mesh deformation solver uses the same mappings as the global solver on the multigrid levels (or is there a reason to use different mappings on the levels? do we transfer the mesh deformation onto the multigrid levels at all? if not, maybe using degree 1 mappings on the multigrid levels is fine? @tjhei can you help here?)
• 3. Since this affects all ASPECT models with free surface we should make sure to have a suitable set of benchmark models and tests to make sure we are not breaking some existing applications. At a minimum all our mesh deformation cookbooks should still work (allken_et_al_2012, continental_extension, crustal_deformation, fastscape_eroding_box, crustal_model_2D/3D, free_surface, free_surface_with_crust). Even before that we should check that relevant benchmarks, in particular the rayleigh_taylor_instability_free_surface, and crameri_et_al benchmarks still give accurate results (ideally better than before, but some are a qualitative code comparisons so hard to do in them).

@tiannh7 I know this is a lot to put on you alone, and we should discuss how to share this work. But since this PR touches one of ASPECTs core algorithms we need to make sure we are not fixing one problem and simultaneously breaking a lot of other models. Can you let us know what you can realistically contribute to the list of points above and we can see where we can help? If there is no chance to make it through all the testing for now, we can also choose to make this improvements optional and defaulting to the old settings, but I agree with the others it would be best to automatically switch everyones models to the better choices.

[gassmoeller]

Two questions here:

• if the model uses higher order velocities but we are in a box geometry we still use a degree 1 mapping? Couldnt we use mesh_deformation_fe.degree to include all possible information we have from the mesh deformation in the mapping? Or do I misunderstand how this works?

• if the model has curved boundaries so far we are using a mapping of degree 4 (in undeformed models). this is on purpose higher than the degree 2 of the velocity FE, because we have seen that spherical geometries improve significantly in accuracy when switching from degree 2 to 4 (I vaguely remember that degree 3 was similar to degree 2, likely this is somehow related to how a spherical geometry would be expanded in polynomials, but I dont have proof and only vague memories of these tests ~10 years ago). Why did you choose mesh_deformation_fe.degree+1 as new degree? And could we use max(4,mesh_deformation_fe.degree+1) instead?

I think my questions show I am a bit uncertain about the relation between the degree of the mapping, and the degree of the vector field that controls the deformation of the geometry. Do they need to be equal? Clearly not, for a default velocity degree of 2 you use mapping degree 3 or 1, depending on the geometry. Does it make sense to use a higher order for the mapping than for the deformation? Probably only if we have additional information from the underlying geometry manifold. Does it make sense to use a lower order for the mapping than for the deformation? I think we are loosing information then, because the lower order mapping cannot represent the higher order deformation.

@tjhei do you have insights here on the optimal degree of the mapping? Or do we need to benchmark?

[gassmoeller]

I am ok either way. The Assert seems slightly too cautious, but using a Stokes degree 1 is quite exotic anyway (since it needs additional stabilization for the finite element to even solve the Stokes equation). So everyone who knows what they are doing and select a velocity degree of 1 likely also knows how to comment out this assert if they are convinced that is what they want.

[gassmoeller]

could you add two more empty lines here?

[gassmoeller]

Nice improvement!

[gassmoeller]

I discovered the need for this in #6284 and did not extensively test or measure it. But: It is done before any particles are created and at this point it is just associating the particle handler with this mapping and triangulation. So at worst we duplicate the initialization of one (or more) empty particle handlers. I wouldnt expect it to be expensive.

[gassmoeller]

So for the level mappings in the GMG solver you always use mesh_deformation_fe.degree, but for the global mapping you use this->get_geometry_model().has_curved_elements() ? (mesh_deformation_fe.degree+1) : 1. I think this may be a problem.

In particular I have found that while the GMG mesh deformation solver is generally faster than the AMG mesh deformation solver, there are some setups where it has serious problems to converge.

E.g. for this parameter file (this is the same you reported Stokes solver issues for, but I got rid of them by increasing the resolution of the model by 1):

geoid_freesurface.prm.txt

I have run the model once with the GMG mesh deformation solver and once with the AMG mesh deformation solver (by commenting out the gmg solver in lines 542-544 in mesh_deformation/interface.cc). The GMG solver is faster as expected, but in timestep 6 and 7 it takes forever to solve, even though it only needs the same number of iterations as the other timesteps. Mesh deformation now takes >90% of the total compute time of this model.

I attach the two log files:
GMG mesh deformation: log.txt

AMG mesh deformation: log.txt

I did not observe the same problems for the free surface cookbook.

I think it is worth looking into the Chebyshev solver settings of the GMG mesh deformation solver if these convergence issues persist, but first we should settle on the optimal mapping degrees (for curved / non curved geometries), then make sure all GMG level mappings use the same degree, and only afterwards check stability/convergence of this test.

[tiannh7]

Hi @gassmoeller,

Thanks a lot for your detailed questions and insights! Here are my thoughts on the points you raised:

I think we need to identify good degrees for the mapping. Intuitively I would think this->get_geometry_model().has_curved_elements() ? std::max(4,mesh_deformation_fe.degree) : mesh_deformation_fe.degree would make sense. However this is different from the current choice, and I dont know how much testing you already did @tiannh7, maybe your choice works better? Why did you choose the numbers you have at the moment?

For my part, I have only tested on some spherical shell examples. In these, I did not observe a significant difference between using degree+1 (i.e., 3) and 4. That said, I agree that using 4 is a more robust choice and likely safer overall.

Then we should make sure the GMG mesh deformation solver uses the same mappings as the global solver on the multigrid levels (or is there a reason to use different mappings on the levels? do we transfer the mesh deformation onto the multigrid levels at all? if not, maybe using degree 1 mappings on the multigrid levels is fine? @tjhei can you help here?)

Yes, I also noticed this issue in the original PR #6284, so I introduced an optional parameter Mesh deformation polynomial degree to give users control over it. However, this wasn’t very transparent to some users. To improve compatibility and simplify things, I switched to using the Stokes velocity polynomial degree as the default, but I overlooked updating this part accordingly.

Since this affects all ASPECT models with free surface we should make sure to have a suitable set of benchmark models and tests to make sure we are not breaking some existing applications. At a minimum all our mesh deformation cookbooks should still work (allken_et_al_2012, continental_extension, crustal_deformation, fastscape_eroding_box, crustal_model_2D/3D, free_surface, free_surface_with_crust). Even before that we should check that relevant benchmarks, in particular the rayleigh_taylor_instability_free_surface, and crameri_et_al benchmarks still give accurate results (ideally better than before, but some are a qualitative code comparisons so hard to do in them).

I fully agree with you on this. I plan to run these tests and benchmarks thoroughly. If my time becomes limited, I will share my results with the team so others can help continue the verification.

Thanks again for your thorough review and helpful suggestions!

[tiannh7]

I have run the model once with the GMG mesh deformation solver and once with the AMG mesh deformation solver (by commenting out the gmg solver in lines 542-544 in mesh_deformation/interface.cc). The GMG solver is faster as expected, but in timestep 6 and 7 it takes forever to solve, even though it only needs the same number of iterations as the other timesteps. Mesh deformation now takes >90% of the total compute time of this model.

Thank you @gassmoeller very much for your detailed feedback and for pointing out the potential performance issues with the GMG solver during later timesteps.

Following your suggestion, I have adjusted the Chebyshev smoother parameters by increasing the smoothing range, polynomial degree, and the number of eigenvalue estimation iterations:

smoother_data[level].smoothing_range = 20.;
smoother_data[level].degree = 7;
smoother_data[level].eig_cg_n_iterations = 15;

While these adjustments brought some improvement, unfortunately, the initialization phase of the smoother still struggles with highly deformed meshes around timesteps 6 and 7, and the solver can still experience significant delays.

I encountered a similar issue earlier in my own 3D spherical domain model, but in my case, the CG iteration count was very high rather than a stall in the smoothing phase. I experimented with different solvers and found that GMRES converged more reliably than CG. This led me to suspect that the problem might be related to the presence of a null space, since for symmetric positive definite Laplace problems, CG should typically converge quickly.

However, the fact that the AMG solver can still work efficiently even under free slip boundary conditions strongly suggests that the physical system itself is well-posed (or at least solvable). This observation directly challenges my initial hypothesis that the null space is the fundamental cause. If the null space were a fundamental issue, AMG should theoretically encounter similar difficulties. Therefore, the problem might be more subtle than a straightforward null space issue.

For this particular case, I also noticed that changing the bottom boundary condition from free slip (set Tangential velocity boundary indicators = inner) to no slip (set Zero velocity boundary indicators = inner) naturally resolved the issue.

Some time ago, I read the thesis by @tjhei’s student Pengfei Jia, Null Space Removal in Finite Element Discretizations, and had two face-to-face discussions with him. However, we did not reach a definitive conclusion.

In short, for 2D/3D spherical models with a free surface on top and free slip at the bottom, the question remains whether null space removal methods like

subsection Nullspace removal
  set Remove nullspace = angular momentum
end

can perfectly eliminate the null space? If there are any defects or corrections needed here, it seems that enforcing a no-slip bottom boundary condition might be a more robust workaround.

Pengfei and I initially planned to reach out to @tjhei for further advice, but due to work commitments, I have delayed contacting him, waiting for a better opportunity. Regardless, these problems are very interesting and the discussions around them have been incredibly helpful for me during my nearly one year working with ASPECT.

I wonder if @tjhei might have any insights or suggestions regarding this issue?

Thank you again for your valuable insights!

[tiannh7]

Hi @gassmoeller,

Sorry for the late reply!

I have recently run tests on the Rayleigh-Taylor instability (RTI) and Crameri case1 (see the attached images) and found that the current improvements have little visible impact on non-curved meshes. This is expected since a linear mapping is already sufficiently accurate for non-curved geometries.

Do you need me to perform any additional tests? Please feel free to raise any questions or concerns.

[gassmoeller]

Hi @tiannh7 thanks for your reply and checking the benchmarks!

The code looks good, I just have two remaining comments on the changelog and test files.

However, I was running some tests with the free_surface cookbook and the fastscape_eroding_box cookbook, and both of them crash with this branch. Here is a screen short of the free_surface cookbook (it is the simpler one), the top panel shows results with the main branch, the bottom shows results with this branch (at the same time step). In the next time step the model with this branch crashes. It looks like either our free surface stabilization needs to be changed for this branch, or there is something else we didnt think about and observe so far.

And it seems the reference tester produces different test results from the ones you included in the PR, could you update the test results with the file from the reference tester jammy-dealii-9.6? Then I can check how large the difference is to the current results.

However, the fact that the AMG solver can still work efficiently even under free slip boundary conditions strongly suggests that the physical system itself is well-posed (or at least solvable). This observation directly challenges my initial hypothesis that the null space is the fundamental cause. If the null space were a fundamental issue, AMG should theoretically encounter similar difficulties. Therefore, the problem might be more subtle than a straightforward null space issue.

I think your idea about the nullspace is a good one, and that AMG works does not necessarily mean much, we extract the constant modes to let the AMG preconditioner handle nullspaces more efficiently. I also noticed that this issue is not directly related to this PR, I had the same issues on a recent development branch without your changes. I will open a separate issue about this, and we dont have to worry about it for this PR.

So I think in order to proceed with this PR we need to figure out what is causing these instabilities I described above and find a fix or workaround. Or we need to leave this optional.

[gassmoeller]

Suggested change

	by ensuring the mesh deformation is consistent with the Stoeks discretization.
	by ensuring the mesh deformation is consistent with the Stokes discretization.

[tiannh7]

Thank you!

[gassmoeller]

We try to avoid including .vtu (or .pvtu, .pvd) files as test output. .vtu files are compressed binary output, so if they change in the future we never know if the change is big or small from the output on Github (without visual inspection in Paraview). It looks like this test relied on the screen-output and statistics files before, and you added the dynamic topography and geoid files (which work nicely, since they are plain text). Please either: 1. remove the visualization output files (.vtu, .pvtu, .pvd, and .visit), or 2. switch the data format of the visualization output to gnuplot (which is plain text), and include the gnuplot output files. Either is ok with me.

[tjhei]

. It looks like either our free surface stabilization needs to be changed for this

I am not 100% sure, but the current stabilization term shouldn't depend on the FE degree.

The free-surface example currently uses theta = 0.5. Would it be worth trying a few different values to see if it makes a difference?

[gassmoeller]

No easy luck, here are versions of the free surface cookbook with default theta (0.5), and maximum stabilization (1.0):

And here is a plot of the maximum topography over time for these two models, and the main branch:

maximum_topography.pdf

From what I gather the stabilization has no effect (it even causes a slightly earlier crash). This could indicate it is not the Stokes velocity that causes the issues (as the stabilization is only applied to the Stokes system, not the mesh deformation system).

I dont know if it is relevant that the crash happens when the maximum topography is reached, at those times the surface is almost not moving.

Alternatively, the issues could also be our explicit time stepping scheme for the mesh deformation (see e.g. this paper that discusses implicit surface deformation schemes instead: https://www.sciencedirect.com/science/article/pii/S0010465515000569)

[tiannh7]

Hi @gassmoeller and @tjhei,

Sorry for the late reply! I just got back home and am currently on the Spring Festival holiday.

Thanks for your detailed testing. It seems the instability arises mainly for Cartesian geometries with free surface and Q2 mesh deformation. While Q2 usually provides better accuracy, it may become unstable under large deformation in this case (which might be caused by the cumulative mesh quality degradation over about 70 steps, rather than due to large deformation in a single time step).

I’ve tried a fallback to Q1 mapping for these cases to improve stability, and visually everything looks fine so far. Do you think this is a reasonable approach, or would you suggest another direction? Or perhaps we could make the mapping order configurable as an option?

unsigned int mapping_degree = mesh_deformation_fe.degree;
...
if (this->get_geometry_model().has_curved_elements())
    mapping_degree = std::max(4u, mapping_degree);
else
    mapping_degree = 1u;

I think your idea about the nullspace is a good one, and that AMG works does not necessarily mean much, we extract the constant modes to let the AMG preconditioner handle nullspaces more efficiently. I also noticed that this issue is not directly related to this PR, I had the same issues on a recent development branch without your changes. I will open a separate issue about this, and we dont have to worry about it for this PR.

Thank you very much for clarifying this! This issue has puzzled me for quite some time. I look forward to the new issue discussion and hope it will lead to a solid resolution.