geodynamics / aspect

A parallel, extensible finite element code to simulate convection in both 2D and 3D models.
https://aspect.geodynamics.org/
Other
227 stars 237 forks source link

Chunk+free surface unwanted deformation #3359

Closed fionaclerc closed 4 years ago

fionaclerc commented 4 years ago

When running the free surface and 2-D 90-degree chunk together, it seems as if the free surface is moving out from the bottom (east) boundary, by ~kms. This is surprising as density/temperature are constant such that no deformation should occur. I'm attaching a plot of the topography at different time steps, as well as a snapshot of paraview showing the pressure field, and the input parameter file. I did play around with different aspect ratios which didn't resolve this, as well as increasing the radial depth which generated instabilities. Thanks for your help ! -Fiona


play_traction6B.prm.txt

pressure_chunk_free_surface

chunk_free_surface_topo.pdf

gassmoeller commented 4 years ago

Hi Fiona, Weird, I do not see anything particular off about this model setup, but I would check the following things:

fionaclerc commented 4 years ago

Hi Rene, Thanks for your reply and suggestions. I did try them:

  • Does this also happen if you use a spherical shell model with the same setup (of course the elements will have a bad aspect ratio, but maybe it still runs)

This leads to small oscillations in the velocity field spherical_shell_velY that correspond to the size of the cells in the direction along the surface (regardless of the inner boundary condition). I realized I get a similar result when using the chunk geometry and imposing a no-slip velocity condition along the inner boundary (instead of free-slip). Putting in a temperature anomaly does override the velocity oscillations; although these oscillations are still visible in the non-adiabatic pressure (which is important as I would like to look at changes in pressure). chunk_NLP

  • Was there a reason you switched to first order temperature elements for this model? Does anything change with 2nd order?

Nothing changes with 2nd order when I was holding the temperature constant; with the temperature anomaly the local velocity field is slightly different but otherwise everything else is similar.

  • Is the nonadiabatic pressure more or less constant this way, or are there oscillations? If so do they affect the small pressure spikes at the top and bottom boundary or does the whole pressure anomaly oscillate between left and right boundary?

There are oscillations between left/right (that decay in magnitude through time) but the non-adiabatic pressure is always negative on the left boundary and positive on the right. The small pressure spikes along the bottom also vary in magnitude accordingly but do not change position.

Is the solution to use the no-slip bottom boundary condition and a very fine mesh ?

Thanks, Fiona

gassmoeller commented 4 years ago

Hi Fiona, I looked into this further and it is a complicated issue (took me the better part of 2 days to figure out the reason, found some unrelated problems with the chunk and still have no perfect solution to the original problem). But please try the workaround I proposed in #3377. There are still the small oscillations that you described for the spherical shell, but I am not sure those are real, and not just a postprocessing artifact. They are pretty small (1e-5 m/yr), and they disappear if you do not interpolate the output solution, i.e. set Interpolate output = false for the visualization postprocessor.

fionaclerc commented 4 years ago

Hi Rene, Thanks for writing the workaround ! The results do look good, although I am having trouble running the same setup with the side boundaries set to the initial lithostatic pressure condition. The sides of the free surface rides up until the model does not converge. Is this a related issue ?

chunk_free_surface_psides.pdf

gassmoeller commented 4 years ago

Is this a related issue ?

I am not sure, but I think not. Can you open another issue with a parameter file?