[WIP] Improve free surface accuracy

[gassmoeller]

This PR is for @tiannh7 to test free surface accuracy.

Do not review yet, this needs some testing and cleanup.
This is based on #4379.

[tiannh7]

Hi Gassmoeller,

I would like to report an error I encountered while running ASPECT 3.1.0-pre with deal.II 9.7.0-pre in a 3D MPI parallel environment (64 processes), in debug mode. The error occurs when resuming from a snapshot, during the solution transfer step (solution_transfer.cc), specifically at line 99 where an exception is thrown due to a data size mismatch:

(data_range.size() % bytes_per_entry == 0) violated

This triggers an abort with a detailed stack trace pointing to dealii::SolutionTransfer and its unpack_callback function. The error indicates that the serialized data size does not match expected sizes during unpacking degrees of freedom values.

Additional remarks:

• The problem does not occur when using the original ASPECT main branch (stable version).
• The error is also not present in 2D chunk geometry models.
• This aligns with what you previously mentioned regarding the stability differences.
• The issue seems related to serialization/deserialization of solution data during mesh refinement or repartitioning in 3D with particles.
• Given the error points to deal.II's SolutionTransfer internals, it might be associated with compatibility or handling of data during solution transfer in the pre-release versions of ASPECT and deal.II.
• I am providing the full debug log for your reference.

Would you have any suggestions or insights on whether this is a known problem, or how to proceed with debugging or workaround?

Thank you very much for your help and time!

particles_error_debug.log

[gassmoeller]

Thanks for the report. I dont have an obvious candidate for what is going wrong, but will think a bit about it. Just to make sure: This checkpoint you are resuming from was already written with the new version of ASPECT and deal.II, correct? Old checkpoints will not be compatible with this version of ASPECT and deal.II (because the field in which we store the mesh displacements is different between this version and the current development version on the main branch).

[tiannh7]

Thank you for your message. Yes, the checkpoint I am resuming from was created with the latest branch of ASPECT and deal.II. I have also just tried running the simulation completely from scratch with this version, but I encountered the same error. Please let me know if you need any further information or logs to help diagnose the issue.

[gassmoeller]

I took a quick look through the code, but didnt see anything obvious that would break with this PR. Could you try to simplify your model setup as much as possible, and then post a simple .prm file that creates the error for you? Then I could try to reproduce the problem.

[tiannh7]

Thank you very much for taking the time to review the code!

I appreciate your suggestion to simplify the model setup. I have prepared a pr-6284.prm.txt file that reproduces the issue I encountered. This parameter file uses ASPECT version 3.1.0-pre (pr-6284, commit 0e7ecaf) together with deal.II 9.7.0-pre (pr-18396, commit 6e3761e52d).

Please let me know if you need any further information or assistance in reproducing the problem. I look forward to your feedback!

[gassmoeller]

Thanks for posting the input file 👍. I will take a look during the next days when time permits. I also found another issue with this PR, which is that models with particles do not work correctly with the changes. Both issues will have to be fixed or at least understood before we can move forward with this PR.

[tiannh7]

Thanks for posting the input file 👍. I will take a look during the next days when time permits. I also found another issue with this PR, which is that models with particles do not work correctly with the changes. Both issues will have to be fixed or at least understood before we can move forward with this PR.

Dear @gassmoeller,

Thank you very much for your attention and time regarding the issues.

I would like to share that I have added a new function in my dealii branch set_mapping_particle_handler to address the particle handler mapping update:

template <int dim, int spacedim>
void ParticleHandler<dim, spacedim>::set_mapping(const Mapping<dim, spacedim> *new_mapping)
{
  mapping = new_mapping;
}

I have not submitted a pull request yet because I have not fully reviewed the dealii contribution guidelines and coding style requirements.

In ASPECT, after the line

sim.mapping.reset(new MappingQEulerian<dim, LinearAlgebra::Vector>...);

I call this new function for each particle manager as follows:

for (auto &pm : sim.particle_managers)
  pm.set_mapping(sim.mapping.get());

Although I am not entirely sure if this is the most appropriate or clean solution, my tests so far indicate that it works correctly without obvious problems.

I sincerely appreciate your ongoing efforts on this, and I look forward to any suggestions or guidance you might have.

[gassmoeller]

Ah, you have thought this through further than I did. Your solution would probably work, but adding new deal.II functions will complicate the transition a bit, because ASPECT needs to stay backward compatible to earlier deal.II for a while. Could you try if using the function ParticleHandler::initialize instead would work as well? That would delete all existing particles, but I hope that at the time we change the mapping we have not yet created the initial particle distribution.

[gassmoeller]

@tiannh7 I uploaded a new version of this PR that seems to restart successfully with the .prm file you provided. Could you give this new version a try? Keep in mind I have not tested this for accuracy or stability, I only made sure it runs. So benchmarking is advised before relying on the results of this changed version.

[tiannh7]

I appreciate your work and will share feedback once testing is complete.

[tiannh7]

@tiannh7 I uploaded a new version of this PR that seems to restart successfully with the .prm file you provided. Could you give this new version a try? Keep in mind I have not tested this for accuracy or stability, I only made sure it runs. So benchmarking is advised before relying on the results of this changed version.

Hi @gassmoeller,

Thank you very much for the updated version of the PR. I have completed the testing and can confirm that the initialization and model restart are working perfectly with the provided prm file.

However, I have not yet performed any accuracy or stability tests.

Also, if I may mention, regarding the visco-elastic problems, I recall from browsing the Use correct stresses in visco-elasto-plastic rheology PR that there seem to be some discrepancies between the particle method and the field method results. This might be something worth keeping in mind during further validation.

Thanks again for your work on this!

[tjhei]

We will need to modify

aspect/source/mesh_deformation/interface.cc

Lines 963 to 966 in ee7942d

	Assert(mesh_deformation_fe.degree == 1, ExcNotImplemented());
	// To be efficient, the operations performed in the matrix-free implementation require
	// knowledge of loop lengths at compile time, which are given by the degree of the finite element.
	const unsigned int mesh_deformation_fe_degree = 1;

to pick the correct degree. I am happy to do the change.

[tiannh7]

Hi Timo,

Thank you very much for your reply and for following up on this.

My preliminary tests have shown that both AMG and GMG exhibit divergence or lack of convergence issues in some non-Cartesian problems. It seems that the entire AMG and GMG workflow might need a thorough review, or this could be related to the EV or SUPG stabilization schemes.

I’ll keep you updated if I make any further progress in my testing.