Initial topography causes mesh holes

[Minerallo]

Hi, As @anne-glerum already mentioned it in the issue #3402, using the initial topography can caused the mesh to have holes in it (see figures). From the tests that I ran, it seems that if the element's size between the elements of the topography and the the one of the internal mesh are too different. The mesh got messed up after re-refining and coarsening some elements. I attached a prm file making the error for a simple case. test_holes.txt

• Initial

• After refinement at timestep 5

NB=In my subduction model, the problem caused locally higher strain rate close to the holes, lower viscosity and higher velocity (my model has diffusion and dislocation). The only way I was able to get rid of the anomalies, was to set a high minimum refinement level (here Level6) to have a resolution of my mesh close to the highest elements resolution of my model (so it becomes over-resolved). A screenshot of the ‘strain rate’ allows to easily notice the anomalies.

Cheers Michael Pons

[gassmoeller]

Hi Michael, thanks for reporting this. This reminds me strongly of problems we had in the past with the spherical geometry. In that case there was no real "hole" in the mesh, it simply appeared that way in the output, because the output vertices were interpolated linearly between neighbour vertices. Just to verify this a bit further, can you test the following:

• Can you repost the second figure if you 'set Interpolate output = false' in the visualization postprocessor?
• Do you use a free surface? If not how does the mesh look like if you do? If yes, how does the mesh look like if you do not? (the reason I am asking is that we use different mappings internally, and that might influence the output)
• You say the strain rate artifacts go away "if I use a high enough resolution", but to me it looks more like the artifacts go away if you use a uniform mesh (i.e. once your minimum refinement is equal to your maximum refinement). Do I see that correctly? Can you check if you see the strain rate artifacts when you have a global uniform mesh of level 4 or 5?

[Minerallo]

Hi Rene, I ran a few tests,

1) I set Interpolate output = false. This model uses free surface.

The anomaly still appears.

• Initial time step

• After rerefinement at timestep =3

2) I made a model where the top is free slip (set Tangential velocity boundary indicators = top)

• Timestep=2(similar to initial timestep), then the model crashed at timestep =3 when it had to re-refine.

• Here is the crash message

*** Timestep 3:  t=6.5265 years

--------------------------------------------------------
An error occurred in line <850> of file </home/ponsm/aspect/source/simulator/core.cc> in function
    void aspect::Simulator<dim>::compute_current_constraints() [with int dim = 2]
The violated condition was: 
    current_constraints.is_consistent_in_parallel( dof_handler.locally_owned_dofs_per_processor(), locally_active_dofs, mpi_communicator, false)
Additional information: 
    Inconsistent Constraints detected!

Stacktrace:
-----------
#0  /home/ponsm/buildaspect/debugfscp2/aspect: aspect::Simulator<2>::compute_current_constraints()
#1  /home/ponsm/buildaspect/debugfscp2/aspect: aspect::Simulator<2>::start_timestep()
#2  /home/ponsm/buildaspect/debugfscp2/aspect: aspect::Simulator<2>::run()
#3  /home/ponsm/buildaspect/debugfscp2/aspect: void run_simulator<2>(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, bool, bool, bool)
#4  /home/ponsm/buildaspect/debugfscp2/aspect: main

3) I set a uniform square mesh of 8x8 km for my model

• Mesh at initial time step

• Mesh after re-refinement, timestep =7

• Viscosity at timestep 7

and strain rate at timestep 7 (sorry for the rainbow colormap)

So it seems that has soon as elements are not uniform anymore after the rerefinement timestep, the "holes appear"

NB : I wanted to use initial topography to avoid the isostatic equilibrium steps, but I found another way by using a time dependent function for velocity boundary. The model first runs a few hundred thousands years without any velocity for the time to reequilibriate and afterward I start pushing.

Cheers, Michael

[gassmoeller]

Hi Michael, I tried reproducing your problem, but was so far unsuccessful. I used the parameter file you attached to the first post, but for me the mesh is coarsened uniformly after timestep 5 (not only in special regions like your figures) and therefore I do not see any holes. I tried with the current master and deal.II 9.1 and the current development version of deal.II. I attach my log file so you can check for any differences to your tests. How many processes did you use? Did you use exactly the prm you attached for the first post? Did you use a different ASPECT or deal.II version?

log.txt

[Minerallo]

Hi Rene, I took back the prm file from the blogpost and ran it again. All the simple case models were run with 5 MPI, I got exactly the same problem. But when I ran it with 1 MPI like you did, there is no problem anymore.

NB : The prm file attached in the first post is the same than the figure but with squared mesh and normally leads to the error. This is a log at the time I ran the expample of the first blog post for a model with squared elements log.txt

So, it seems that it is a matter of how many number of processes I am using. Using the same file tried : -1MPI works -2MPI works -3MPI Problem -4MPI works -5MPI Problem -6MPI Problem for only 1 timestep after refinement -7MPI Problem -8MPI works -9MPI Problem -10MPI Problem

this is what I am using now to run the MPI tests -- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion. -- . version 2.2.0-pre (old_version, 145bba1) -- . using deal.II 9.1.1 -- .with 32 bit indices and vectorization level 1 (128 bits) -- . using Trilinos 12.10.1 -- . using p4est 2.2.0 -- . running in DEBUG mode -- . running with 3 MPI processes the log : log.txt

Bests, Michael

[tjhei]

We just discussed during the ASPECT biweekly meeting. This might be due to a missing call to Triangulation::communicate_locally_moved_vertices() https://www.dealii.org/developer/doxygen/deal.II/classparallel_1_1distributed_1_1Triangulation.html#a247598f1323a9f847832e60d6c840469

That said, I don't think this covers the case when we refine and repartition completely, but only to make sure the ghost layer is consistent. Is there a way to redo the operation for initial topography after every refinement step? If not, I propose making initial topography a mesh deformation plugin instead #3399. This would make this problem disappear.

[ricitron]

I noticed this problem as well in 2d and 3d spherical chunk geometry with initial topography. The hole is hard to see in the image below but it is just to the right of the xyz orientation axis where the refinement gets coarser.

Is this a real hole or an issue with the plotting? I am using set Interpolate output = false.

The hole occurs where the change in slope of the surface topography intersects a coarser cell at depth.

Here is a case with more extreme topography showing more holes:

[bangerth]

I think that @tjhei 's hunch is probably right, but it's also possible that the mesh movement uses a way of moving nodes individually without respecting that the positions of hanging nodes cannot be chosen independently of their surrounding node locations.

[ricitron]

This should be resolved now with merged pull request 3864

[gassmoeller]

Since this seems resolved I will close the issue for now. Feel free to reopen if the problem persists.